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1.  SUMMARY 


In  this  project  (FA9453-10-C-0263)  we  build  on  our  earlier  results  (DE-FC52- 
06NA27324,  Dreger  et  al.,  2008;  Ford  et  al.,  2010;  Ford  et  al.,  2009a, b;  Ford  et  al.,  2008) 
to  investigate  whether  source-type  populations  (explosion,  collapse,  earthquakes)  that 
were  found  to  separate  for  the  Western  US  remain  separated  for  other  regions  of  the 
world.  We  have  expanded  our  database  of  studied  events  to  include  natural  and  man¬ 
made  seismicity  from  the  Korean  Peninsula,  Kazakhstan,  and  China.  For  Kazakhstan  and 
China  we  have  performed  waveform  modeling  to  develop  calibrated  ID  velocity  models 
to  improve  regional  distance  moment  tensor  capabilities  in  the  region,  and  have  refined 
an  approach  that  combines  the  use  of  long-period  regional  distance  waveforms  and  short- 
period  P-wave  first-motions.  This  combined  approach  is  found  to  improve  discrimination 
capabilities  under  very  sparse  seismic  monitoring  conditions.  In  addition,  we  have 
compiled  source-type  distributions  for  natural  seismicity  for  California  for  use  as  a 
reference  in  the  development  of  a  statistical  discrimination  approach.  We  have  developed 
a  continuous  scanning  method  that  can  be  applied  to  targeted  regions  for  autonomous 
event  detection,  location,  source  estimation  and  discrimination. 

2.  INTRODUCTION 

Ford  et  al.  (2009a)  calculated  seismic  moment  tensors  for  17  nuclear  test  explosions, 
12  earthquakes,  and  3  collapses  in  the  vicinity  of  the  Nevada  Test  Site  in  the  Western  US. 
They  found  that  the  relative  amount  of  isotropic  and  deviatoric  moment  provided  a  good 
discriminant  between  the  explosions,  earthquakes  and  collapses.  Importantly,  the 
populations  of  explosions  and  cavity  collapses  separate  dramatically  whereas  in 
traditional  mb:MS  and  regional  phase  P/S  ratio  discriminants  this  separation  is  less  clear 
(e.g.  have  inconsistent  results  (Walter  et  al.,  2007)). 

In  this  report  we  present  the  analysis  of  several  nuclear  and  industrial  explosions  to 
extend  the  database  of  explosion  moment  tensor  solutions.  In  particular  we  have 
developed  solutions  for  the  1988  Joint  Verification  Experiment  event  at  the 
Semipalantinsk  test  site  in  Kazakhstan,  two  Chinese  tests  at  the  Lop  Nor  test  site  that 
occurred  in  1995  and  1996,  and  an  industrial  explosion.  We  have  applied  a  standard 
seismic  moment  tensor  method  (Minson  and  Dreger,  2008),  the  network  sensitivity 
solution  (NSS)  (Ford  et  al.,  2010),  and  the  combination  of  long-period  complete 
waveforms  and  P-wave  first  motions  NSS  method  (Ford  et  al.,  2012),  which  is  described 
in  section  3.  The  work  has  resulted  in  three  peer-reviewed  publications  (Ford  et  al.,  2012; 
Chiang  et  al.,  2014;  Nayak  and  Dreger,  2014). 
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3.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1.  Event  Discrimination  using  Regional  Moment  Tensors  with 
Teleseismic-P  Constraints 

Introduction 

The  elements  of  the  seismic  moment  tensor  can  be  used  to  derive  the  source  type. 
Determination  of  the  source  type  and  discrimination  between  explosion,  collapse,  and 
earthquake  sources  can  aid  in  the  understanding  of  deformation  in  a  region  and  in  nuclear 
explosion  monitoring.  Ford  et  al.  (2009b)  calculated  a  seismic  moment  tensor  for  the 
announced  nuclear  test  of  the  Democratic  People’s  Republic  of  Korea  (DPRK)  on  25 
May  (Memorial  Day)  2009  using  regionally  recorded  (<1500  km)  intermediate  period 
(10-50  s)  waveforms.  The  source  type  derived  from  the  seismic  moment  tensor  was 
dominantly  explosive.  An  earthquake  source  was  highly  unlikely  due  to  the  poor  fit  to  the 
data  for  any  double-couple  (DC),  or  earthquake-like,  source.  However,  another  source 
type  with  very  little  explosive  component  fit  the  data  almost  equally  as  well  as  an 
explosion. 

This  source  is  described  as  a  compensated  linear- vector  dipole  (CLVD)  with  a 
vertical  axis  in  compression  (VCLVD-P).  The  similarity  in  data  fit  between  the 
dominantly  VCLVD-P  source  and  dominantly  explosive  source  presents  a  problem  in 
discrimination  between  these  two  source  types.  The  reason  for  this  VCLVD-P-explosion 
ambiguity  is  that  the  regional  distance  intermediate  period  recordings  of  a  shallow  event 
are  most  sensitive  to  the  isotropic  tensile  radiation  pattern  along  the  equator  of  the  focal 
sphere,  which  produces  Rayleigh  waves  with  uniform  source  phase,  and  which  does  not 
generate  Love  waves.  The  pure  explosion  and  VCLVD-P  radiation  patterns  differ 
significantly  for  a  small  takeoff  angle  where  waves  leave  the  source  nearly  vertical.  The 
incorporation  of  additional  data  sensitive  to  this  vertical  region  of  the  focal  sphere  would 
constrain  the  VCLVD-P-explosion  tradeoff.  We  propose  to  use  teleseismic-P  recordings 
to  constrain  the  moment  tensor  derived  source  type  because  the  steep  takeoff  angle  of 
teleseismic-P  waves  is  sensitive  to  the  lower  hemisphere  of  the  focal  sphere,  where  they 
are  in  dilation  for  a  VCLVD-P  and  in  compression  for  an  explosion. 

The  use  of  only  teleseismic  data  to  determine  the  source  type  has  more  serious  trade¬ 
offs  in  terms  of  explosion  monitoring,  where  there  is  little  discrimination  between 
shallow  reverse  mechanism  earthquakes,  explosions,  and  VCLVD  in  tension.  However, 
in  the  absence  of  regional  data,  relative  amplitudes  of  teleseismic  body  waves  offer  some 
ability  to  determine  source  type,  and  this  was  the  approach  of  Clark  and  Pearce  (1988), 
Pearce  et  al.  (1988),  and  Rogers  and  Pearce  (1992).  The  combined  use  of  regional  and 
teleseismic  data  to  constrain  the  source  type  was  done  qualitatively  by  Bowers  (1997)  to 
discriminate  between  a  collapse,  explosion,  and  an  earthquake  in  a  South  African  mine 
and  by  Bowers  and  Walter  (2002)  to  discriminate  between  a  collapse  and  an  explosion 
for  mine  events  in  Volkershausen,  Germany,  and  Wyoming,  United  States.  In  those 
studies  the  teleseismic-P  waveform  was  shown  to  be  inconsistent  with  sources 
determined  via  regional  waveform  fitting.  A  similar  approach  is  used  here,  but  with  a 
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more  quantitative  component  afforded  by  the  use  of  a  source-type  plot  (Hudson  et  al., 
1989)  with  application  to  the  2009  Memorial  Day  explosion  in  the  DPRK. 

Figure  1  shows  the  distribution  of  regional  and  far-regional  distance  stations  that  were 
used  in  the  analysis,  and  the  global  locations  of  three  arrays  that  were  used  to  obtain  P- 
wave  polarities.  Seismic  moment  tensor  solutions  were  obtained  using  the  Network 
Sensitivity  Solution  (NSS)  method  developed  by  Ford  et  al.  (2010).  In  Figure  2  the 
waveforms  for  several  mechanisms  sampled  in  the  NSS  space  are  compared.  It  is  notable 
that  while  there  is  essentially  the  same  level  of  fit  between  the  full  moment  tensor 
solution  and  the  deviatoric  moment  tensor  solution  there  is  a  marked  difference  in  fit 
level  with  double-couple  solutions.  The  NSS  source-type  space  (Figure  2c)  shows  the 
theoretical  tradeoff  between  a  explosion  source  and  a  CLVD  with  a  vertical  major  vector 
dipole  in  compression  (a  negative  CLVD).  The  NSS  represents  the  best  fit  surface  of  a 
uniform  distribution  of  moment  tensor  solutions  in  source-type  space.  It  is  important  to 
note  that  the  construction  of  the  NSS  involves  the  forward  testing  of  upwards  of  30 
million  moment  tensor  solutions,  and  that  there  are  many  mechanisms  of  different 
orientations  (orientation  of  eigenvectors,  strike/dip/rake,  and  azimuth  and  plunge  of 
CLVD  major  vector  dipole). 

The  three  dots  in  each  of  the  focal  mechanisms  show  the  locations  of  the  P-wave 
observations  derived  from  array  analysis  at  three  global  distance  arrays.  The  first  motions 
are  all  up.  It  is  evident  that  the  first  motions  eliminate  some  of  the  double-couple 
solutions,  however  the  regional  long-period  waveforms  clearly  eliminate  the  possibility 
of  such  solutions  for  this  event.  The  first  motions  are  also  seen  to  be  inconsistent  with  the 
down  first  motion  region  of  the  deviatoric  solution. 

Starting  with  the  NSS  shown  in  Figure  2c  we  tested  solutions  against  the  teleseismic 
first  motion  observations,  where  if  a  particular  moment  tensor  solution  violates  the  first 
motion  observations  the  goodness  of  fit  parameter  for  that  solution  is  set  to  zero.  Figure 
2d  shows  the  results  of  this  test.  The  deviatoric  line  in  the  source  type  diagram  is 
eliminated  from  possible  solutions,  and  the  best  fitting  region  is  tightly  defined  around 
the  full  moment  tensor  inversion  result.  The  use  of  teleseismic  first  motions  eliminates 
the  theoretical  tradeoff  between  pure  explosion  and  negative  CLVD  sources. 

Figure  2c  and  2d  demonstrate  that  the  incorporation  of  teleseismic-P  information  can 
eliminate  the  theoretical  tradeoff  between  explosion  and  -CLVD  sources  due  to 
similarity  in  surface  wave  radiation  patterns,  and  greatly  improve  the  characterization  of 
the  explosive  nature  of  the  event.  This  work  was  published  in  the  Bulletin  of  the 
Seismological  Society  of  America  (Ford  et  al.,  2012). 

In  section  3.2  we  refine  this  approach  to  utilize  regional  P-wave  first  motion  data,  as 
well  as  a  modified  goodness  of  fit  criterion. 
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Figure  1.  Map  of  the  Yellow  Sea  and  Korean  Peninsula  with  the  2009  Memorial  Day 
explosion  (star)  and  the  regional  stations  used  in  the  analysis  (triangles).  Inset:  azimuthal 
equidistant  map  of  teleseismic  arrays  used  in  the  analysis  (triangles),  where  the  origin  is 
the  location  of  the  explosion. 


4 


Approved  for  public  release;  distribution  is  unlimited. 


(a) 


Models,  and  regional  waveforms 


Full  Dev  DCJ 

•  00 


DC2  Explosion 


£ 

3 


< 

x 


(b)  Teleseismic-P  wayelrain 


(c)  Regional  NSS 


EiploSLtHl 


Implosion 


(d)  Regional  NSS  with  teleseismic-P  constraint 
Explosion 


Figure  2.  Teleseismic  constraint  to  regional  source-type  determination,  (a)  Three- 
component  (tangential,  T;  radial,  R;  and  vertical,  V)  intermediate-period  (10-5 Os) 
displacements  (labeled  Data)  from  regional  stations  shown  in  Figure  2  and  synthetic 
waveforms  at  those  stations  predicted  by  the  sources  shown  with  a  lower-hemisphere 
equal-area  focal  sphere.  Amplitudes  are  normalized  to  the  maximum  of  the  three 
components  at  each  station.  The  gray  bar  near  the  bottom  is  50  s  long.  The  gray  scale  of 
the  compressional  component  of  the  focal  sphere  matches  the  gray  scale  of  the  predicted 
waveform  and  is  also  noted  for  MDJ  waveforms.  The  T  and  P  axes  are  plotted  on  the 
focal  sphere  of  the  best-fit  full  moment  tensor  (labled  Full).  The  dots  show  the  teleseismic 
P-wave  polarity  (all  up)  for  the  three  teleseismic  arrays,  (b)  Array  beams  of  high- 
frequency  (0.8  to  4.5  Hz)  teleseismic-P  wave  trains  (labeled  Data).  As  in  (a)  the  gray 
scale  relates  predicted  solutions  to  specific  source  types,  (c)  Contours  of  fit  (variance 
reduction),  VR,  equation  1)  to  the  regional  data  for  a  uniform  distribution  of  sources  on  a 
source-type  plot,  (d)  The  NSS  shown  in  (c)  where  VR  is  now  zero  if  the  teleseismic-P  and 
forward  predicted  data  are  uncorrelated,  as  is  the  case  for  the  deviatoric  source. 
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3.2.  Analysis  of  the  1988  Soviet  JVE  and  1990’s  Chinese  Lop  Nor  Tests 

Introduction 

The  use  of  regional  distance  long-period,  complete  waveform  data  to  determine 
the  seismic  moment  tensor  and  discriminate  the  source-type  of  earthquakes,  underground 
cavity  collapse  and  nuclear  explosions  has  been  demonstrated  for  events  in  the  western 
United  States  (Dreger  et  al.,  2008;  Ford  et  al.,  2008;  Ford  et  al.,  2009a),  and  for  the 
recent  2006  and  2009  North  Korean  nuclear  tests  (Ford  et  al.,  2009b;  Ford  et  al.,  2010). 
In  these  studies  populations  of  earthquakes,  underground  cavity  collapses  and  nuclear 
explosions  are  found  to  separate  when  considered  on  a  Hudson  et  al.  (1989)  source-type 
diagram.  The  source  type  plot  simplifies  the  moment  tensor  into  two  parameters  that 
depend  on  the  eigenvalues  of  the  moment  tensor.  These  parameters  T  and  k  describe  the 
deviation  from  a  pure-double-couple  in  terms  of  non-volumetric  (compensated  linear 
vector  dipole,  CLVD)  and  volumetric  components.  Ford  et  al.  (2010)  utilized  the  Hudson 
et  al.  (1989)  source-type  representation  to  develop  a  network  sensitivity  solution  (NSS), 
which  tests  on  the  order  of  100  million  moment  tensor  solutions  uniformly  distributed  in 
source-type  space  to  determine  the  best  fitting  solution,  the  uncertainty  in  the  solution, 
and  the  capabilities  of  the  method  given  the  topology  of  a  recording  station  network.  The 
regional  distance  moment  tensor  inversion,  coupled  with  NSS  analysis,  and  the 
characterization  of  sensitivities  and  uncertainties  due  to  random  errors  and  systematic 
velocity  model  errors  enable  the  discrimination  of  source-type  in  relatively  sparse 
regional  distance  monitoring. 

We  investigate  the  14  September  1988  US-Soviet  Joint  Verification  Experiment 
(JVE)  nuclear  test  at  the  Semipalatinsk  test  site  in  Eastern  Kazakhstan,  and  two  nuclear 
explosions  conducted  less  than  ten  years  later  at  the  Chinese  Lop  Nor  test  site  (Table  1). 
These  events  were  very  sparsely  recorded  by  stations  located  within  1600  km,  and  in 
each  case  only  3  or  4  stations  were  available  in  the  regional  distance  range  (Figure  3  a)  for 
moment  tensor  analysis.  Following  the  results  of  Ford  et  al.  (2009b)  we  incorporated 
first-motion  data  from  the  regional  stations,  as  well  as  teleseismic  stations  to  provide 
additional  constraint  in  the  NSS  analysis.  The  results  show  that  unique  discrimination  of 
these  events  is  possible  under  these  extremely  sparse  monitoring  conditions  when  long- 
period  regional  waveforms  and  P-wave  first-motion  polarities  are  combined. 
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Table  1.  Event  Hypocenter  Information 


1988  JVE 

1995  Lop  Nor 

1996  Lop  Nor 

1999  Lop  Nor 
EQ 

Origin  Time 

03:59:57.30 

04:06:00.20 

02:56:00.06 

03:51:05.42 

(UTC) 

Latitude, 

49.882, 

41.553, 

41.5804, 

41.674,  88.463 

longitude 

78.882 

88.7496* 

88.6893* 

Depth  (km) 

<1.0 

0.7* 

0.5* 

17 

Mbf 

6.1 

6.1 

5.9 

5.9 

MsJ 

4.8 

5.0 

4.3 

5.3 

Mw 

5.2 

5 .2-5.4 

5.2 

5.3 

F 

0.97  to  1.15 

3.79  to  4.20 

-12.75  to  -20.44 

K 

2.59  to  2.73 

2.80  to  2.88 

2.81  to  2.90 

T 

-0.87  to  -0.91 

-0.82  to  -0.83 

-0.62  to  -0.66 

0.1424 

K 

0.57  to  0.59 

0.56  to  0.57 

0.36  to  0.39 

0.0540 

*  Waldhauser  et  a 

[.  (2004) 

f  Priestly  et  al.  (1990);  Yang  et  al.  (2003) 

%  ISC  catalog  (Explosions);  NEIC  catalog  (Earthquake) 

Data  Processing  and  Methods 

Data  and  instrument  response  for  the  US-Soviet  Joint  Verification  Experiment 
(JVE)  and  Lop  Nor  nuclear  tests  were  downloaded  from  the  Incorporated  Research 
Institutions  for  Seismology  (IRIS,  see  Data  and  Resources  section)  that  consists  of  a 
collection  of  stations  from  the  regional  broadband  seismic  network  in  China,  the  Global 
Seismic  Network  (GSN),  and  other  temporary  networks  (Figure  3a).  In  addition  to  data 
from  IRIS  for  the  1988  Soviet  JVE,  we  also  have  data  from  temporary  deployments  of 
short  period  instruments  located  at  near-regional  distances  <  260  km  (Priestley  et  al., 
1990)  and  the  Borovoye  station  (BRVK)  at  a  distance  of  690  km.  The  data  was 
instrument  corrected,  integrated  to  displacement,  rotated  to  radial  and  transverse 
components,  and  filtered  between  12.5  to  80  seconds  with  a  Butterworth  filter  depending 
on  instrument  type  and  signal-to-noise  levels.  Table  2  lists  the  stations  used  for  each 
event,  and  the  frequency  passband  that  was  employed.  The  processed  data  were  then 
inverted  using  the  time  domain  full  waveform  moment  tensor  inversion  of  Minson  and 
Dreger  (2008).  For  JVE  we  used  a  source  depth  of  1  km  that  gives  the  highest  goodness 
of  fit  between  data  and  synthetics,  and  for  the  two  Lop  Nor  events  we  used  the  source 
depths  from  Waldhauser  et  al.  (2004).  We  then  inverted  the  data  with  Green’s  functions 
computed  for  a  range  of  source  depths  to  determine  the  source  depth  sensitivity. 
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Figure  3.  a)  Event  locations  (star)  and  seismic  stations  (triangle  and  square)  used  in  the 
moment  tensor  analysis.  The  two  Lop  Nor  events  are  located  very  closely  hence  the 
overlapping  stars.  Triangles  represent  the  stations  used  in  the  analysis  of  the  15  May 
1995  Lop  Nor  explosion,  8  June  1996  Lop  Nor  explosions,  and  squares  are  the  stations 
used  in  the  14  September  1988  Soviet  Joint  Verification  Experiment.  Focal  mechanisms 
of  local  earthquakes  used  in  the  velocity  model  calibration  are  also  plotted.  On  the  right 
is  the  solution  from  this  study  and  on  the  left  is  the  Harvard  GCMT  solution. 

Figure  3.  b)  Velocity  models  used  in  the  analysis.  BAY  from  Walter  and  Ammon  (1993); 
MODI  from  broadband  waveform  modeling  of  local  earthquakes  (Fig.  3a);  MOD2  is  a 
ID  model  simplified  from  Sun  et  al.  (2010)  ’s  surface  wave  tomography;  NE  is  a  simple 
layer  over  halfspace  model. 

The  seismic  moment  tensor  consists  of  nine  force  couples  that  represent  the 
equivalent  body  forces  for  seismic  sources  of  different  geometries  (Jost  and  Herrmann, 
1989),  that  due  to  conservation  of  angular  momentum  reduce  to  six  independent  couples 
and  dipoles.  The  data  are  represented  by  the  convolution  of  Green’s  functions  for  a  given 
Earth  model,  source  terms  and  the  moment  tensor  elements.  The  individual  moment 
tensor  elements  are  obtained  using  a  generalized  least  square  inversion  and  the  goodness 
of  fit  between  the  data  and  synthetics  is  measured  by  the  variance  reduction  (VR)  given 
by: 


VR  = 


1- 


IW-*)2 

z* 


TOO 


(1) 


where  d  represents  the  data  and  s  represents  the  synthetic  waveforms. 

Because  the  data  are  linear  combinations  of  the  Green’s  functions  weighted  by 
their  associated  moment  tensor  elements,  we  need  a  well-calibrated  velocity  model  in 
order  to  estimate  robust  seismic  source  parameters.  For  each  test  site  we  utilized  several 
published  seismic  velocity  models,  and  for  some  paths  we  additionally  calibrated  models 
by  modeling  earthquake  records.  The  inversion  method  also  allows  for  small  time  shifts 
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between  the  data  and  Green’s  functions  to  compensate  for  errors  in  origin  time,  location, 
and  in  velocity  structure. 

In  addition,  we  use  two  alternate  methods  in  order  to  fully  characterize  the  source 
solution  space.  One  utilizes  a  grid  search  to  find  the  best-fitting  double-couple  (DC), 
pure-isotropic/explosion  (ISO),  or  DC  +  ISO  source  mechanism.  The  other,  described 
below,  utilizes  a  moment  tensor  grid  search  in  order  to  assess  solution  uniqueness  and 
resolution  of  source-type. 

Although  the  seismic  moment  tensor  inversion  gives  a  unique  moment  tensor 
solution,  the  decomposition  of  the  moment  tensor  solution  is  non-unique.  Therefore  we 
implemented  the  Hudson  et  al.  (1989)  source-type  representation  that  circumvents  the 
need  to  decide  on  a  particular  moment  tensor  decomposition  scheme.  The  source-type 
diagram  has  two  key  parameters  T=-2s  and  k  on  the  x-  and  y-axis,  respectively,  given  by 
the  equations: 


-m.  •  „ 

s= — 71,  m3  +  m2+m]  =  0, 


m. 


m. 


>m2> 


m. 


(2) 


K  = 


\MJ  + 


m. 


M  =(m3+m2+m1) 


5  iso 


(3) 


m'i  and  m'3  are  the  deviatoric  principal  moment  associated  with  the  minimum  and 
maximum  principal  compressive  stress  axes,  and  Miso  is  the  isotropic  moment.  Equation 
(2)  measures  the  deviation  from  a  pure  shear  dislocation  and  equation  (3)  describes  the 
volume  change.  In  this  convention,  s  is  0  for  a  pure  DC  source  and  ±0.5  for  a  pure  CLVD 
source,  and  k  is  ±1  for  a  spherical  explosion  and  implosion,  respectively.  Understanding 
the  relative  contributions  of  the  different  moment  tensor  elements  provides  insights  into 
the  complex  source  processes  of  explosions  as  well  as  other  seismic  events.  This 
representation  of  the  seismic  source  has  been  shown  to  result  in  separate  populations  for 
explosions,  underground  cavity  collapse  and  earthquakes  (Ford  et  al.,  2009a;  Ford  et  al., 
2009b;  Dreger  et  al.,  2008;  Ford  et  al.,  2008)  enabling  discrimination  capability. 
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Table  2.  Even-Station  Combinations  &  Passbands 


Event 

Station 

Passband  (s) 

Distance  (km) 

Azimuth  (deg) 

1988  JVE 

BAY 

12.5-20 

255 

257 

KKL 

12.5-20 

256 

295 

BRVK 

12.5-30 

690 

304 

WMQ 

20-50 

953 

132 

1995  Lop  Nor 

MAK 

30-50 

796 

319 

AAK 

30-50 

1184 

281 

TLY 

30-50 

1597 

40 

1996  Lop  Nor 

MAK 

30-50 

794 

319 

AAK 

30-80 

1184 

280 

TTB02 

30-50 

1410 

246 

1999  Lop  Nor  EQ 

WMQ 

30-50 

246 

346 

PDG 

30-50 

760 

287 

MAKZ 

30-50 

770 

320 

TLG 

30-50 

940 

284 

TLY 

30-50 

1602 

41 

To  assess  the  confidence  of  the  moment  tensor  solution,  we  implemented  the 
Network  Sensitivity  Solution  (NSS)  technique  developed  by  Ford  et  al.  (2010).  The 
technique  presents  the  level  of  fit  between  actual  data  and  the  different  theoretical 
solutions  described  by  the  source-type  diagram  for  a  given  station  configuration,  Earth 
model,  and  frequency  band.  From  the  NSS  of  a  given  event  we  can  determine  whether  or 
not  the  best  fitting  full  moment  tensor  solution  from  the  inversion  is  well-resolved  to 
make  useful  interpretations  about  the  source.  We  included  regional  and/or  teleseismic  P- 
wave  first  motions  in  addition  to  waveform  data  in  the  NSS  analysis  (Ford  et  al.,  2012)  to 
better  constrain  the  moment  tensor  solution.  To  include  P-wave  polarities  as  additional 
constraints,  we  take  the  suite  of  synthetic  moment  tensor  solutions  from  the  previous 
waveform  NSS,  compute  their  P-wave  polarities,  and  compare  them  to  the  observed  P- 
wave  polarities.  We  assign  -1  for  downward  motion  and  +1  for  upward  motion.  The  VR 
is  calculated  as: 


VR  = 


T.POIL 


TOO 


(4) 


Then  we  calculate  the  combined  waveform  and  first  motion  VR  as: 


VR  =  (.v  ■  .v  VRfi3Mkmi )  ■  1 00  (5) 

where  sVRreg  and  sVRfm  are  normalized  by  the  maximum  regional  waveform  and  first 
motion  VR,  respectively.  The  take-off  angles  are  calculated  using  iaspei-tau  (Snoke, 
2009)  and  the  iasp91  reference  Earth  model.  Incorporating  the  first  motion  data  proves  to 
be  a  powerful  tool  eliminating  solution  non-uniqueness,  and  in  assessing  the  confidence 
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of  a  solution  under  sparse  station  coverage  monitoring  situations.  We  find  that  additional 
polarity  constraints  assist  by  uniquely  discriminating  the  events  as  predominantly 
explosive. 

This  powerful  result  is  due  to  the  fact  that  the  tradeoff  between  a  pure  explosion 
and  negative  CLVD  (Ford  et  al.,  2012)  only  occurs  for  a  negative  CLVD  in  which  the 
compressional  major  vector  dipole  is  vertically  oriented.  For  these  two  sources  and 
combinations  of  the  two,  the  surface  wave  radiation  patterns  are  identical  for  Rayleigh 
waves,  and  are  theoretically  null  for  Love  waves.  This  vertical  CLVD  mechanism 
predicts  dilational  first  motions  in  the  center  of  the  focal  sphere.  For  explosions,  the 
incorporation  of  P-wave  first  motion  data,  particularly  teleseismic  observations  with  low 
takeoff  angle  proves  to  be  an  effective  test  against  the  negative  CLVD  source  type. 
Naturally  challenges  remain  for  cases  with  low  signal-to-noise  in  P-waves,  or  cases  in 
which  ffee-surface  reflections,  or  interactions  with  nearby  velocity  structure  result  in 
reversed  dilational  polarities.  Furthermore,  the  absence  of  good  P-wave  polarities  from 
teleseismic  arrays  or  stations  may  require  the  use  of  regional  seismic  waveform  data 
alone.  Nevertheless,  as  will  be  shown  in  the  following,  even  in  such  regional  distance, 
sparse  coverage  situations,  the  combination  of  low-frequency  full  waveform  fits  with 
regional  distance  P-wave  first-motions,  greatly  enhances  the  discriminatory  power  of  the 
method. 

14  September  1988  Soviet  Joint  Verification  Experiment  (JVE) 

Our  analysis  of  the  Soviet  JVE  event  included  three  stations  in  Eastern 
Kazakhstan  and  one  station  in  northwestern  China.  Stations  BAY  and  KKL  are  filtered 
between  12.5  to  20  seconds,  station  BRVK  is  filtered  between  12.5  to  30  seconds,  and 
station  WMQ  is  filtered  between  20  to  50  seconds.  The  use  of  station-specific  filters 
allows  tailoring  of  the  approach  for  path  specific  signal-to-noise  levels,  and  with  respect 
to  the  suitability  of  the  velocity  model  used  to  compute  the  Green’s  functions.  We  used  a 
well-calibrated  velocity  model  BAY  (Walter  and  Ammon,  1993)  for  the  Eastern 
Kazakhstan  stations  and  tested  two  models  for  the  path  to  the  Chinese  station  WMQ.  The 
two  models  used  for  station  WMQ  are  from  broadband  waveform  modeling  of  two  local 
earthquakes  in  Kazakhstan  and  the  Lop  Nor  nuclear  test  site,  called  MODI  (Fig.  3b),  and 
another  ID  model  simplified  from  surface  wave  tomography  (Sun  et  al.,  2010),  called 
MOD2.  Figure  3b  shows  all  the  velocity  models  used  in  this  paper  to  calculate  the 
Green’s  functions  and  generate  the  synthetic  seismograms.  We  have  both  regional  P- 
wave  first  motions  including  the  four  stations  used  in  the  waveform  inversion  (Walter 
and  Patton,  1990;  Priestley  et  al.,  1990),  and  teleseismic  P-wave  first  motions  from  the 
Adirondack  array  in  New  Hampshire,  reported  by  Battis  and  Cipar  (1991),  and 
Gauribidanur  array  in  India.  The  best  fitting  source  mechanism  from  full  moment  tensor 
inversion  of  regional  waveforms  consists  of  a  predominantly  explosive  component  and  a 
moment  magnitude  (Mw)  between  5.21  and  5.25,  depending  on  the  velocity  model  used. 
The  calculated  moment  used  the  convention  for  total  moment  described  by  Bowers  and 
Hudson  (1999).  Mw  falls  close  to  the  surface  wave  magnitude  (Ms)  in  the  catalog 
because  the  frequency  range  we  used  in  the  inversion  is  dominated  by  long  period  surface 
waves.  Respective  VRs  for  a  full,  deviatoric,  DC,  DC+ISO,  and  ISO  mechanism  are  84- 
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85%,  79-80%,  75-76%,  81-82%,  and  77-78%.  Variations  in  VR  resulted  from  different 
velocity  models  (MODI  and  MOD2)  used  for  WMQ.  The  goodness  of  fit  between  data 
and  synthetics  are  similar  for  the  five  different  moment  tensor  decompositions  shown  in 
Figure  4,  due  to  the  dominance  of  surface  waves,  and  presence  of  Love  waves,  but  a  pure 
DC  solution  has  the  poorest  VR  as  compared  to  the  other  four  decompositions  that  are 
either  purely  isotropic  or  included  an  isotropic  component. 
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Figure  4.  1988  Soviet  JVE  Source  Model  and  Waveform  Comparisons.  The  top  trace  is 
the  observed  waveform  data  and  the  rest  are  synthetic  waveforms  according  to  their 
corresponding  source  mechanisms  (The  order  from  top  to  bottom  is:  data,  full, 
deviatoric,  double-couple  (DC),  DC  +  isotropic  (ISO),  and  ISO).  Station  name,  source- 
receiver  distance,  azimuth  and  maximum  amplitude  are  shown.  Associated  variance 
reduction  (VR)  for  each  source  type  is  plotted  next  to  their  focal  mechanisms.  Data  shows 
small  Love  waves  on  the  tangential  component,  characteristic  of  an  explosion.  We  show 
results  using  Green ’s  functions  calculated  from  MOD2  for  WMQ. 

We  computed  the  NSS  using  three  different  combinations  of  data  sets:  (1)  from 
regional  waveforms,  (2)  combined  waveform,  regional  and  teleseismic  P-wave  first 
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motions,  and  (3)  combined  waveform  and  four  regional  P-wave  first  motions.  We  also 
compare  the  NSS  using  different  Earth  models  MODI  (Fig.  5a-c)  and  MOD1+MOD2 
(Fig.  5d-e)  to  compute  the  Green’s  functions.  The  regional  waveforms  only  NSS  solution 
(Fig.  5a, d)  shows  a  similar  trend  compare  to  other  nuclear  explosions,  such  as  the  NTS 
events  and  the  two  DPRK  tests  (Ford  et  al.,  2010),  with  the  best-fitting  full  moment 
tensor  solution  plotting  near  the  theoretical  opening  crack  mechanism  and  there  being  a 
trend  in  the  best  fitting  region  extending  to  the  negative  CLVD.  The  shaded  contour 
regions  correspond  to  different  scaled  variance  reduction  (sVR)  in  which  the  sVR  in 
Figure  5  is  scaled  to  the  moment  tensor  solution  in  the  NSS  that  has  the  maximum  VR. 
Moment  tensor  solutions  fitting  >  90%  sVR  covers  the  upper  right  half  of  the  Hudson 
source  type  plot,  whereas  solutions  fitting  >  98%  sVR  wraps  around  a  small  region 
around  the  theoretical  opening  crack  and  includes  the  best  fitting  mechanism  from  the 
full  moment  tensor  inversion.  In  the  case  of  using  just  the  waveform  data,  source 
mechanisms  without  a  significant  explosive  component  can  fit  the  observed  data  just  as 
well  as  a  predominantly  explosive  mechanism  (Fig.  5a, d).  However,  when  regional  and 
teleseismic  P-wave  (Fig.  6)  first  motions  are  included  in  the  computation  of  the  NSS  a 
solution  that  is  predominately  explosive  is  obtained  (Fig.  5b, e).  The  NSS  results  show 
significant  improvement  in  discrimination  capabilities  when  we  included  additional 
constraints  from  P-wave  first  motions,  especially  for  moment  tensor  solutions  fitting 
better  than  sVR  of  90%.  Figure  5c, f  suggests  including  just  the  regional  P-wave  polarity 
measurements  from  the  same  four  stations  used  in  the  waveform  inversion  can  greatly 
increase  discrimination  capabilities,  which  is  good  since  good  teleseismic  data  may  not 
always  be  available.  The  additional  constraints  from  P-wave  first  motions  eliminate  the 
common  ISO-CLYD  tradeoff  as  well  as  the  mechanisms  that  do  not  agree  with  both  the 
higher  frequency  polarity  data  and  the  long  period  waveforms. 


Figure  5.  Network  Sensitivity  Solutions  (NSS)  for  the  1988  Soviet  JVE.  Circles  are  source 
mechanisms  from  Figure  4;  crosses  are  the  theoretical  mechanisms;  shaded  regions  are  full 
moment  tensor  solutions  contoured  according  to  their  scaled  variance  reduction  (sVR).  Best 
fitting  full  and  deviatoric  mechanisms  are  also  plotted.  Black  crosses,  open  circles  and  circles 
plotted  on  top  are  P-wave  up  first  motions,  T-  and  P-axes.  a-c)  NSS  using  MODI,  d-f)  NSS  using 
MODI  and  MOD2. 
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15  May  1995  Lop  Nor  Shaft  Explosion 

In  contrast  with  the  1988  Soviet  JVE,  the  Lop  Nor  events  studied  here  had  only 
three  available  regional  stations  with  adequate  signal-to-noise  in  the  intermediate-  to 
long-periods  employed  by  our  method,  and  therefore  this  case  presents  a  very  sparse 
monitoring  scenario.  We  used  broadband  waveform  data  filtered  between  30  to  50 
seconds  period  from  three  stations  in  Eastern  Kazakhstan,  Eastern  Kyrgyzstan  and 
Central  Siberia.  For  the  furthest  station  TLY  in  Siberia  we  used  a  simple  layer  over 
halfspace  velocity  model,  NE.  This  model  was  obtained  using  forward  and  inverse 
modeling  of  a  local  earthquake  in  Lop  Nor,  the  same  earthquake  used  to  obtain  MODI. 
We  used  the  iasp91  crustal  velocity  as  a  starting  model,  and  from  our  modeling  result  we 
observe  that  because  of  the  distance  and  long  periods  the  waves  are  not  sensitive  to  the 
finer  details  of  the  velocity  model,  and  additional  layers  to  the  halfspace  model  do  not 
improve  the  waveform  modeling  results  significantly.  For  station  AAK  in  Kyrgyzstan  we 
used  MODI,  and  for  station  MAK  in  Kazakhstan  we  tested  both  MODI  and  MOD2.  The 
original  instrument  response  file  for  MAK  from  IRIS  resulted  in  anonymously  low 
waveform  amplitudes.  After  comparing  the  coda  envelope  functions  with  nearby  stations 
we  concluded  the  original  response  file  could  not  be  used  because  the  amplitude  of  the 
coda  envelope  function  for  MAK  is  significantly  shifted  from  all  other  nearby  stations, 
indicative  of  a  problematic  instrument  response.  Rautian  and  Khalturin  (1978)  observed 
that  coda  envelope  functions  decay  stably  over  time,  and  amplitudes  only  vary  from  event 
to  event.  Instead  we  used  a  response  file  from  the  previous  date  for  the  analysis.  The  first 
motion  picks  for  the  1995  Lop  Nor  event  are  all  from  regional  stations  obtained  through 
IRIS,  including  the  three  stations  used  in  the  moment  tensor  inversion  (Fig.  6). 
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Figure  6.  P-wave  first  motions  for  the  earthquake  and  two  explosions.  The  waveforms  are 
bandpass  filtered  between  0.2  to  3  Hz,  except  BRVK  for  the  1988  JVE,  which  is  filtered 
between  0.2-1. 5  Hz  to  avoid  exceeding  the  Nyquist  frequency.  P-wave  polarities  at  KSU 
and  the  Adirondack  array  data  not  shown  here  (used  in  the  JVE  analysis)  are  from 
Figure  3  in  Priestley  et  al.  (1990)  and  Figure  12  in  Battis  and  Cipar  (1991),  respectively. 
All  waveforms  are  in  velocity  except  for  the  1988  JVE  event,  which  is  in  displacement. 
Map  shows  the  distribution  of  stations  and  array  beams  for  the  three  events. 
Gauribidanur  array  beam  (GBB)  is  filtered  between  0.8-5  Hz. 

The  best  fitting  full  moment  tensor  is  a  predominantly  explosive  mechanism  and 
has  comparable  VR  to  the  deviatoric,  DC  and  DC+ISO  mechanisms,  but  a  pure  ISO 
source  does  not  fit  well.  Respective  VRs  for  a  full,  deviatoric,  DC,  DC+ISO  and  ISO 
mechanisms  are  81-86%,  76-78%,  73-74%,  74-78%,  and  10%,  depending  on  the  velocity 
model  used  for  MAK.  A  pure  explosive  source  cannot  generate  Love  waves  to  fit  the 
large  observed  Love  wave  amplitudes  on  the  tangential  component,  and  therefore  the 
pure  explosion  model  has  the  anonymously  low  VR  of  10%  (Fig.  7a).  This  is  of  concern 
since  natural  tectonic  earthquakes  represented  by  a  DC  actually  fit  better  than  the  pure 
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explosion,  however  in  similarity  to  the  JVE,  a  pure  DC  solution  has  the  poorest  VR 
compared  to  other  decompositions  that  include  an  isotropic  component  (excluding  a  pure 
ISO  source).  Difficulty  fitting  the  relatively  small  Rayleigh  waves  on  radial  and  vertical 
components  at  AAK  is  likely  caused  by  Love  wave  multipathing  based  on  the  particle 
motions.  Surface  waves  on  the  radial  and  vertical  components  exhibit  horizontally 
polarized  Love  wave  particle  motion.  The  best  fitting  full  moment  tensor  solution  is 
located  near  the  opening  crack  mechanism  on  a  Hudson  source  type  plot.  This  result  is 
consistent  with  the  JVE  event  and  previously  studied  nuclear  explosions  (Lord  et  al., 
2010).  The  waveform  only  NSS  results  show  a  wide  range  of  possible  sources  fitting  > 
90%  of  the  best  fitting  moment  tensor  solution  (Fig.  7b, d),  which  is  largely  the  result  of 
the  large  Love  wave  amplitudes  and  the  sparse  station  coverage.  However,  if  we  use  both 
waveform  data  and  P  wave  polarities  observed  at  regional  distances  we  see  the  combined 
analysis  significantly  reduces  the  distribution  of  solutions  with  high  sVR  (>  90%)  and 
uniquely  discriminates  the  event  as  consistent  with  other  nuclear  explosions  and 
inconsistent  with  earthquakes  and  collapses  (Fig.  7c, e).  Figure  7c, e  shows  that  the 
double-couple  and  deviatoric  mechanisms  that  can  fit  the  surface  wave  data  well  fail  to  fit 
the  first  motions  at  the  sites.  Velocity  model  variations  also  show  some  variation  in  the 
NSS  as  expected;  however,  the  distribution  of  sources  fitting  >  98%  of  the  best  fitting 
solution  is  similar  for  both  velocity  models  illustrating  that  the  approach  is  capable  of 
source  type  discrimination  even  in  cases  where  the  earth  structure  is  not  as  well 
constrained. 
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Figure  7.  a)  1995  Lop  Nor  Source  Model  and  Waveform  Comparisons.  The  top  trace  is 
the  observed  waveform  data  and  synthetic  waveforms  according  to  their  corresponding 
source  mechanisms.  Associated  variance  reduction  (VR)  for  each  source  type  is  plotted 
next  to  their  focal  mechanisms.  Data  show  big  Love  waves  on  the  tangential  components. 
We  show  results  using  Green ’s  functions  calculated  from  MODI  for  MAK. 

Figure  7.  b-e)  Network  Sensitivity  Solutions  (NSS)  for  the  1995  Lop  Nor  Explosion. 
Circles  are  best  fitting  full,  deviatoric,  pure  DC,  ISO+DC  and  pure  ISO  source 
mechanisms;  crosses  are  the  theoretical  mechanisms;  shaded  regions  are  full  moment 
tensor  solutions  contoured  according  to  their  scaled  variance  reduction  (sVR).  Crosses, 
open  circles  and  filled  circles  plotted  on  top  are  P-wave  up  first  motions,  T-  and  P-axes. 
Figure  7.  b-c)  NSS  using  MODI;  Figure  7.  d-e)  NSS  using  MODI  and  MOD2. 
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8  June  1996  Lop  Nor  Shaft  Explosion 

For  the  second  Lop  Nor  event  we  used  regional  waveform  data  from  three 
broadband  stations  MAK  in  Kazakhstan,  AAK  in  Kyrgyzstan,  and  TTB02  in  northeastern 
Pakistan,  and  displacement  data  and  synthetics  are  filtered  between  30  to  50  seconds,  30 
to  80  seconds  and  30  to  50  seconds,  respectively.  For  AAK  and  TTB02  we  used  MODI 
to  calculate  the  Green’s  functions,  and  for  MAK  we  used  both  MODI  and  MOD2  since 
its  location  falls  within  the  surface  tomography  study  of  Sun  et  al.  (2010).  We  used  P- 
wave  first  motions  at  regional  distances  and  one  station  at  teleseismic  distances  (Fig.  6). 
All  stations  are  obtained  from  IRIS.  The  qualities  of  the  picks  are  good  and  most  show 
impulsive  upward  P-wave  first  motions  on  the  vertical  components.  The  1996  Lop  Nor 
event  has  a  slightly  smaller  mb  of  5.9  reported  in  the  catalog  compared  to  the  1995  Lop 
Nor  event,  which  has  an  mb  of  6.1.  The  Mw  from  the  regional  moment  tensor  inversion  is 
around  5.2  depending  on  the  velocity  models  used  (Table  1).  Similar  to  that  of  the  1995 
event  we  observe  significant  Love  waves  on  the  tangential  component.  Because  of  the 
strong  Love  waves  the  goodness  of  fit  between  data  and  synthetic  for  a  pure  ISO  solution 
could  not  fit  the  data  well  (<10%  VR)  and  the  best  fitting  mechanism  is  an  implosion 
instead  of  an  explosion.  To  further  explore  how  a  pure  isotropic  solution  fits  the  data,  we 
searched  for  the  best  fitting  isotropic  mechanism  using  the  vertical  displacement  only.  A 
pure  implosion  mechanism  fits  station  TTB02  and  AAK  at  44%  VR  and  51%  VR, 
respectively  with  a  seismic  moment  of  1.7xl022  dyne-cm,  but  does  not  fit  station  MAK, 
in  fact  the  synthetic  is  phased-shifted  by  ±71  from  the  data.  Significant  contribution  from 
non-isotropic  sources  resulted  in  poor  fits  to  the  data  when  we  assume  a  purely  isotropic 
source.  It  is  possible  that  the  tectonic  release  for  this  event  is  large  and  likely  caused 
Rayleigh  wave  reversals  (Toksoz  and  Kehrer,  1972;  Ekstrom  and  Richards,  1994)  thus 
the  pure  ISO  mechanism  is  more  implosive  rather  than  explosive  (Fig.  8a).  Rayleigh 
wave  reversals  due  to  large  tectonic  release  have  been  observed  in  other  nuclear 
explosions  as  well,  such  as  the  1998  Indian  test  (Walter  and  Rodgers,  1999). 
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Figure  8.  a)  1996  Lop  Nor  Source  Model  and  Waveform  Comparisons.  The  top  trace  is 
the  observed  waveform  data  and  synthetic  waveforms  according  to  their  corresponding 
source  mechanisms.  Associated  variance  reduction  (VR)  for  each  source  type  is  plotted 
next  to  their  focal  mechanisms.  Data  shows  big  Love  waves  on  the  tangential  component. 
We  show  results  using  Green ’s  functions  calculated  from  MODI  for  MAK. 

Figure  8.  b-e)  Network  Sensitivity  Solutions  (NSS)  for  the  1996  Lop  Nor  Explosion. 
Circles  are  best  fitting  full,  deviatoric,  pure  DC,  ISO+DC  and  pure  ISO  source 
mechanisms;  crosses  are  the  theoretical  mechanisms;  shaded  regions  are  full  moment 
tensor  solutions  contoured  according  to  their  scaled  variance  reduction  (sVR).  Crosses, 
open  circles  and  filled  circles  plotted  on  top  are  P-wave  up  first  motions,  T-  and  P-axes. 
Figure  8.  b-c)  NSS  using  MODI;  Figure  8.  d-e)  NSS  using  MODI  and  MOD2. 
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The  full  and  deviatoric  solutions  for  the  1996  Lop  Nor  event  have  comparable  VR 
at  78%  but  are  not  significantly  higher  than  the  best  fitting  pure  DC  and  ISO+DC 
mechanisms  that  have  VR  at  73%  (Fig.  8a).  However,  a  vertical  dip-slip  mechanism  is  an 
uncommon  shallow  crustal  earthquake.  The  majority  of  near- vertical  or  sub-horizontal 
dip-slip  mechanisms  in  the  GCMT  catalog  are  associated  with  subduction  and  spreading 
centers.  We  expect  the  waveform  only  NSS  result  cannot  uniquely  discriminate  the  1996 
event  as  a  predominantly  explosive  source  due  to  the  strong  Love  waves  and  sparse 
station  coverage.  This  is  indeed  the  case  and  the  waveform  only  NSS  shows  many 
different  moment  tensor  solutions  have  high  sVR  >  90%,  and  although  the  NSS  contours 
with  sVR  >  98%  includes  mostly  explosion-like  mechanisms  it  also  extends  downward 
and  crosses  into  the  horizontal  deviatoric  axis  that  is  not  observed  in  the  1996  event  (Fig. 
8b, d).  This  behavior  is  likely  due  to  Rayleigh  wave  reversal.  However,  after  we 
incorporated  regional  P-wave  first  motions  the  combined  NSS  results  now  show  similar 
trends  as  observed  in  the  1995  Lop  Nor  test  and  the  1988  Soviet  JVE,  though  the 
difference  is  contours  showing  solutions  with  sVR  >  90%  are  more  extensive  and  cross 
slightly  over  to  the  horizontal  deviatoric  line  (Fig.  8c, e).  Although  the  1996  Lop  Nor 
combined  waveform  and  first  motion  NSS  does  not  give  a  unique  discrimination,  it 
identifies  the  source  as  non-DC.  Unlike  earthquakes,  the  distribution  of  moment  tensor 
solutions  is  not  situated  around  the  pure  DC  mechanism  but  shifted  along  the  vertical 
volumetric  axis  and  towards  an  opening  dipole.  The  best  fitting  full  moment  tensor  for 
the  1996  Lop  Nor  explosion  lies  close  to  the  opening  dipole,  whereas  for  the  Lop  Nor 
earthquake  the  solution  lies  close  to  the  DC.  Although  the  1996  Lop  Nor  event  is  not 
uniquely  discriminated  as  a  predominantly  explosive  source,  NSS  shows  the  event  is 
unusual  and  unlike  a  typical  earthquake  that  needs  further  analysis  to  fully  characterize 
its  source  process.  The  earthquake  shown  in  Figure  9  is  one  of  earthquakes  used  to  obtain 
MODI,  located  close  to  the  Lop  Nor  test  site  and  at  17  km  depth.  Five  stations  were 
available  in  Eastern  Kazakhstan  and  central  Siberia  for  the  inversion  (Fig.  9a)  and  NSS 
(Fig.  9b-c)  but  the  data  are  nosier  compare  to  the  explosion  data. 
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Figure  9.  Waveform  fits  and  Network  Sensitivity  Solutions  (NSS)  for  the  30  January 
1999  Earthquake  near  the  Lop  Nor  Test  Site.  The  distribution  of  solutions  resembles  the 
Little  Skull  Mountain  aftershock  (Ford  et  al.,  2010)  but  the  center  of  the  bullseye  is 
shifted  to  the  right  towards  the  negative  CLVD.  The  station  coverage  is  slightly  better 
than  the  Lop  Nor  explosions  with  additional  stations  to  the  west,  giving  better  azimuthal 
coverage. 


Sensitivity  to  Station  Geometry  and  Source  Depth 

Manmade  explosions  are  conducted  at  shallow  depths.  For  the  1988  Soviet  JVE 
we  used  a  source  depth  of  1.0  km  and  for  the  two  Lop  Nor  shaft  explosions  we  used 
initial  source  depths  reported  in  Waldhauser  et  al.  (2004),  which  were  estimated  using  the 
relationships  between  mb,  yield,  and  source  depth.  To  test  our  method’s  sensitivity  to 
source  depth  we  performed  regional  waveform  moment  tensor  inversion  using  a  suite  of 
Green’s  functions  calculated  at  source  depths  ranging  from  0.2  km  to  16  km.  Unlike 
tectonic  earthquakes,  explosions  have  much  shallower  source  depth,  generally  <  1  km. 
Therefore,  source  depth  can  be  a  very  useful  discriminant  for  earthquakes  and  explosions. 

Plotting  all  of  the  best  fitting  full  moment  tensor  solutions  with  additional 
constraints  from  P-wave  first  motions  at  different  depths  on  a  Hudson  source-type  plot 
(Fig.  10),  we  see  the  sVR  gets  much  worse  as  source  depth  increases.  The  1988  JVE  and 
1995  Lop  Nor  event  solutions  with  high  sVR  and  shallow  depths  (<  1  km)  are  centered  at 
the  theoretical  opening  crack  mechanisms.  In  comparison,  source  depth  for  the  1996  Lop 
Nor  event  is  less  constrained.  Higher  uncertainty  from  the  Green’s  functions  (path 
between  the  event  and  station  TTB02)  and  possibly  greater  tectonic  release  for  this  event 
may  have  resulted  in  less  sensitivity  to  source  depth.  We  did  not  specifically  model  the 
path  between  the  event  and  station  TTB02,  hence  greater  errors  may  be  introduced  into 
the  calculation  of  the  Green’s  functions.  Although  the  1996  Lop  Nor  event  source  depth 
is  not  as  well  constrained  as  the  other  two  events  in  this  paper,  most  of  the  high  VR 
solutions  are  shallow  and  the  mechanisms  are  close  to  an  opening  linear  vector  dipole. 
Solutions  at  greater  depths  are  not  explosive  but  also  have  lower  VR. 
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In  addition  to  depth  sensitivity,  we  also  performed  Jackknife  tests  for  the  three 
events  to  discern  the  inversions’  sensitivity  to  station  geometry.  The  Jackknife  tests 
reveal  we  need  at  least  three  stations  to  have  confidence  in  the  moment  tensor  solution. 
Solutions  obtained  using  less  than  three  stations  are  not  stable,  and  depending  on  station 
geometry  can  give  you  incorrect  source  mechanism.  Generally  moment  tensors  computed 
using  stations  closer  to  the  source  (-300-700  km,  as  in  the  case  of  JVE)  resolve  more 
isotropic  component  and  can  better  constrained  the  NSS,  regardless  of  the  total  number  of 
stations  used  in  the  inversion. 


Figure  10.  Source  depth  sensitivity  for  the  three  events  analyzed.  Best- fitting  full  moment 
tensors  for  the  1988  JVE  (diamond),  1995  Lop  Nor  (square),  and  1996  Lop  Nor  (circle) 
based  on  regional  waveform  inversion  and  P-wave  first  motion  constraints.  Size  of  the 
symbol  increases  with  decreasing  source  depth,  filled  symbols  have  sVR<  90%,  and  open 
symbols  have  sVR  >  90%. 


Effects  of  Tectonic  Release  and  Source  Medium  Damage 

The  presence  of  shear  waves  in  seismic  recordings  of  nuclear  explosions  indicates 
the  explosion  process  is  complex  and  involves  other  source  processes  such  as  interactions 
with  the  ffee-surface,  the  effects  of  shockwave  and  relaxation  of  tectonic  strain  (Patton 
and  Taylor,  2011;  Toksoz  et  al.,  1965;  Toksoz  and  Kehrer,  1972;  Burger  et  al.,  1986; 

Day  et  al.,  1987).  Recent  studies  (Ben-Zion  and  Ampuero,  2009;  Patton  and  Taylor, 
2011)  have  suggested  that  significant  contributions  to  the  radiated  seismic  wavefield  can 
arise  from  source  medium  damage.  Material  damage  associated  with  explosions  affects 
the  Rayleigh  wave  radiation  pattern  (Patton  and  Taylor,  2008),  therefore  damage  may 
have  significant  implications  to  source-type  discrimination.  Our  combined  waveform  and 
first  motion  analysis  of  the  Lop  Nor  tests  shows  the  regional  moment  tensor  is  a 
promising  method  for  source-discrimination  even  when  tectonic  release  is  large.  The  best 
fitting  solutions  located  in  the  region  between  a  pure  explosion  and  tensile  crack  on  the 
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source  type  diagram  is  comparable  to  solutions  for  the  DPRK  (Ford  et  al.,  2009b;  Ford  et 
al.,  2010)  and  NTS  explosions  (Ford  et  al.,  2009a),  and  is  consistent  with  the  nature  of 
tensile  damage  above  the  shot  point  due  to  spall  and  other  free  surface  interactions  as 
proposed  by  Patton  and  Taylor  (2008).  To  further  assess  the  capabilities  of  regional 
moment  tensors  for  events  with  tectonic  release  or  source  damage,  we  conducted  a  series 
of  synthetic  tests  to  examine  these  effects. 

Two  types  of  composite  sources  were  tested:  (1)  a  ISO  +  DC  source  and  (2)  a  ISO 
+  CLVD  source.  The  ISO  +  DC  source  examines  the  effects  of  tectonic  release  and 
consists  of  an  explosion  and  a  45-degree  dipping  reverse  fault  mechanism.  The  ISO  + 
CLVD  source  examines  the  effects  of  source  medium  damage  and  consists  of  an 
explosion  and  a  CLVD  mechanism  with  a  vertically  oriented  major  vector  dipole  (Patton 
and  Taylor,  2008).  To  describe  the  relative  seismic  moments  between  the  different  source 
types,  we  used  the  relationships  proposed  by  Toksoz  et  al.  (1965)  for  tectonic  release  and 
Patton  and  Taylor  (2008,  2011)  for  material  damage  from  a  deviatoric  source.  The  index 
F  measures  the  ratio  between  isotropic  moment  (Miso)  and  double-couple  moment  (Mix:), 


F  =  1.5Mdc 


M iso 

The  index  K  measures  the  relative  strengths  of  the  moment  tensor  elements, 

2Mv 


K  = 


'  xz 
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and  when  tectonic  release  occurs  only  in  the  horizontal  plane  as  in  our  second  test  case, 
2  {Miso  T  Mclvd  ) 

iso  CLVD 


We  used  the  1988  Soviet  JVE  station  geometry  and  Green’s  functions  from  our 
moment  tensor  analysis  to  compute  the  synthetic  data,  and  added  real  noise  collected 
from  the  stations  used  in  the  Soviet  JVE  analysis  or  from  other  nearby  IRIS  stations  when 
pre-event  seismic  recordings  are  not  available.  For  the  two  composite  sources  (ISO+DC 
and  ISO+CLVD)  we  tested  different  F-factors  and  K-factors,  and  for  each  factor  we 
looked  at  different  SNRs  at  5%,  10%,  15%  and  20%  of  the  maximum  vertical  Rayleigh 
wave  amplitude.  We  performed  the  synthetic  testing  using  long  period  waveforms  only, 
and  another  one  combining  the  waveforms  and  four  P-wave  first  motion  measurements 
recorded  at  the  same  seismic  stations.  For  the  remaining  part  of  this  section  we  focused 
on  results  from  the  tests  where  SNR  is  20%  of  the  maximum  vertical  Rayleigh  wave 
amplitude.  This  is  the  most  extreme  case  of  SNR  and  generally  the  real  data  we  use  have 
SNR  levels  that  are  much  lower  than  20%.  The  distribution  of  individual  moment  tensor 
solutions  >  98%  sVR  in  the  Hudson  source-type  diagram  are  more  compacted  for  lower 
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SNR  levels,  but  the  contours  outlining  all  possible  points  are  very  similar  for  all  SNR 
levels.  This  observation  applies  to  both  composite  source  cases  we  tested. 

The  ISO+DC  composite  source  case  shows  the  waveform  only  NSS  becomes 
more  constrained  with  increasing  contribution  from  tectonic  release  (Fig.  1 1).  Best-fitting 
full  moment  tensor  solutions  center  around  the  opening  crack  mechanism,  which  is 
consistent  with  solutions  obtained  from  real  explosions.  Although  this  seems 
counterintuitive,  the  NSS  gives  a  predominantly  explosive  source  because  the 
contribution  from  a  reverse  fault  mechanism  helps  eliminate  the  well-known  ISO-CLYD 
tradeoff  between  explosion  and  -CLVD.  The  reverse  fault  radiation  pattern  is  more 
similar  to  a  +CLVD  mechanism.  Hence,  the  two  competing  mechanisms  eliminate  the 
ISO-CLVD  tradeoff.  Best-fitting  full  moment  tensor  solutions  of  different  degree  of 
tectonic  release  are  all  near  the  opening  crack  mechanism,  which  is  consistent  with 
solutions  obtained  from  real  explosions.  The  waveforms  only  NSS  and  combined 
waveform  and  first  motion  NSS  show  very  similar  distribution  for  sVR  >  98%.  We  see 
no  significant  improvement  by  adding  first  motions  because  the  waveforms  only  NSS 
already  uniquely  identifies  the  source  as  predominantly  explosive.  However,  first  motions 
do  help  constrain  the  NSS  for  solutions  with  sVR  >  90%,  due  to  the  fact  that  the  P-wave 
radiated  from  the  explosion  is  primary  and  precedes  motions  due  to  secondary  block 
faulting  (delay  time  between  explosion  and  secondary  faulting  is  1  second),  or  tensile 
damage  above  the  shot  point.  By  using  both  long  period  waveforms  and  P-wave  polarities 
we  can  capture  the  volumetric  signature  of  the  event  when  tectonic  release  has  significant 
contribution  to  the  explosion  source  processes.  In  the  event  that  we  cannot  uniquely 
discriminate  the  event  as  explosive,  the  deviation  from  deviatoric  mechanisms  flags  the 
event  as  uncommon  and  warrants  further  analysis  to  understand  its  source  processes. 

To  look  at  the  effects  of  medium  source  damage,  we  adopted  the  model  proposed 
by  Patton  and  Taylor  (2008)  where  the  cause  of  damage  is  represented  by  a  vertically 
oriented  CLVD  source.  The  deformation  observed  at  the  source  with  this  model  is 
extensional  along  the  vertical  axis  and  horizontal  compressions  around  the  explosion 
(+CLVD).  As  defined  by  Patton  and  Taylor  (2008),  a  K-factor  of  1  means  contribution 
from  CLVD  vanishes  to  zero,  K  >  1  means  extension  along  the  vertical  axis,  and  K  <  1 
means  compression  along  the  vertical  axis.  Our  synthetic  tests  with  just  waveform  data 
show  K  >  1  (implying  greater  material  damage)  reduces  the  commonly  observed  ISO- 
CLVD  tradeoff  (Fig.  12).  This  is  the  result  of  competing  mechanisms  between  the  - 
CLVD  and  +CLVD,  increasing  contribution  from  a  +CLVD  source  mechanism  affects 
the  Rayleigh  wave  radiation  pattern  thus  enhancing  the  explosive  characteristics  of  the 
source.  In  contrast,  K  <  1  where  compressional  deformation  is  along  the  vertical  axis  the 
ISO-CLVD  tradeoff  persists.  In  some  cases  when  K  is  high  (>  2.5)  the  NSS  identifies  the 
event  as  explosive  even  without  additional  constrains  from  P-wave  first  motions.  In  other 
cases  where  the  event  is  not  uniquely  identified  as  explosive,  additional  P-wave  first 
motions  improve  the  solution  by  shifting  the  98%  sVR  contours  above  the  deviatoric 
CLVD  mechanism  and  reducing  the  region  of  possible  source  mechanisms. 
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Waveforms  Only 


iso  +  DC 


Figure  11.  Network  Sensitivity  Solutions  (NSS)  for  an  ISO+DC  composite  source. 
Plotted  are  the  98%  scaled  variance  reduction  (sVR)  contours  with  increasing  F-factor 
(therefore  increasing  contribution  from  the  DC  component),  and  dashed  contours  are 
from  a  pure  reverse  fault  (DC)  mechanism.  The  circles  are  the  best  fitting  full  moment 
tensor  solution  from  waveform  inversion.  Here  we  compare  how  the  NSS  changes  when 
additional  P-wave  first  motions  are  included. 
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Waveforms  Only  iso  +  CLVD 


Waveforms  and  First  Motions 


Figure  12.  Network  Sensitivity  Solutions  (NSS)  for  an  ISO+CLVD  composite  source. 
Plotted  are  the  98%  scaled  variance  reduction  (sVR)  contours  with  increasing  K-factor 
(therefore  increasing  contribution  from  the  CLVD  component).  The  circles  are  the  best 
fitting  full  moment  tensor  solution  from  waveform  inversion.  Here  we  compare  how  the 
NSS  changes  when  additional  P-wave  first  motions  are  included. 

Conclusions 

We  have  performed  seismic  moment  tensor  inversions  for  the  1988  Soviet  JVE 
test  and  two  Lop  Nor  nuclear  tests.  These  cases  represent  sparse  monitoring  conditions 
and/or  uncertainty  in  velocity  structure.  In  each  case  we  have  shown  that  the  use  of  long- 
period  waveform  data  comprised  mostly  of  regional  surface  waves  result  in  solutions 
with  large  isotropic  components  that  are  consistent  with  solutions  for  other  studied 
nuclear  tests  (Ford  et  al.,  2009a;  Ford  et  al.,  2009b;  Ford  et  al.,  2010).  Using  only 
regional  waveforms,  the  distribution  of  solutions  on  the  source  type  diagram  of  Hudson  et 
al.  (1989)  do  not  cleanly  discriminate  the  event  either  because  of  the  known  explosion 
negative  CLVD  tradeoff  (case  of  the  JVE  event)  or  due  to  large  observed  Love  waves 
(cases  of  the  two  Lop  Nor  tests).  In  each  case,  however,  the  inclusion  of  regional  P-wave 
polarities,  and  ideally  observations  from  teleseismic  arrays  when  available,  reduces  the 


26 

Approved  for  public  release;  distribution  is  unlimited. 


area  of  solutions  that  provide  a  good  level  of  fit  to  the  data,  providing  good  separation 
from  double-couple  solutions  and  solutions  on  the  deviatoric  line. 

The  8  June  1996  Lop  Nor  test  indicates  that  Rayleigh  waves  may  be  reversed  at 
some  of  the  regional  stations  due  to  a  large  tectonic  release.  The  large  tectonic  release  is 
possibly  due  to  shock  driven  block  faulting  or  tensile  damage  (Patton  and  Taylor,  2008). 
To  further  investigate  this  we  carried  out  a  series  of  synthetic  tests  using  the  JVE  station 
geometry  to  evaluate  the  regional  moment  tensor  method  capabilities  for  different  levels 
of  tectonic  release  measured  by  both  the  Toksoz  F  index  and  the  Patton  K  index.  These 
tests  show  that  the  combination  of  long-period  regional  waveform  data  and  regional 
distance  P-wave  first  motions  are  able  to  resolve  the  anomalous  volumetric  nature  of 
compound  explosion  tensile  damage  and  block  faulting  events. 

This  work  has  been  published  in  the  Bulletin  of  the  Seismological  Society  of  America 
(Chiang  et  al.,  2014). 
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3.3.  Analysis  of  the  Effects  of  Vanishing  Traction  on  Seismic  Moment 
Tensor  Recovery  for  Shallow  Explosions 

Introduction 

For  the  nuclear  explosion  source-type  identification  problem  the  uncertainty  in  a 
solution  is  as  important  as  the  best  fitting  parameters.  A  potential  issue  for  shallow 
seismic  sources  that  are  effectively  at  the  free-surface  between  the  ground  and  air  is  that 
the  vanishing  traction  at  the  free-surface  can  cause  the  associated  vertical  dip-slip  (DS) 
Green’s  functions  to  have  vanishing  amplitudes  (Julian  et  ah,  1998),  which  in  turn  can 
result  in  the  indeterminacy  of  the  Mxz  and  Myz  components  of  the  moment  tensor  and 
therefore  bias  in  the  moment  tensor  solution.  The  effects  of  the  free-surface  on  the 
stability  of  the  moment  tensor  becomes  important  as  we  continue  to  investigate  and 
improve  the  capabilities  of  regional  full  waveform  moment  tensor  inversion  for  source- 
type  identification  and  discrimination.  It  is  important  to  understand  its  effects  for 
discriminating  shallow  explosive  sources  for  nuclear  monitoring,  but  could  also  be 
important  in  natural  systems  that  have  shallow  seismicity  such  as  volcanic  environments 
and  geothermal  systems,  and  other  manmade  shallow  seismicity  related  to  anthropogenic 
activities  such  as  hydraulic  fracturing  and  mining. 

To  tackle  the  potential  issues  that  could  arise  in  the  moment  tensor  analysis  of 
shallow  sources,  we  performed  a  series  of  synthetic  tests  to  understand  the  effects  of  free- 
surface  vanishing  traction  on  the  total  seismic  moment,  isotropic  seismic  moment  and  the 
source  mechanism.  We  evaluate  the  sensitivity  of  the  moment  tensor  solutions  as  a 
function  of  source  depth  and  velocity  model.  Based  on  what  we  learn  from  the  synthetic 
studies,  we  applied  the  moment  tensor  method  to  the  HUMMING  ALBATROSS  quarry 
blast  events,  which  is  an  excellent  dataset  in  terms  of  understanding  the  effects  of  free- 
surface  vanishing  traction  in  real  seismic  recordings.  These  small  chemical  explosions  are 
approximately  10  m  deep  and  are  recorded  at  up  to  several  km  distances.  Therefore  the 
data  represent  a  rather  severe  source-station  geometry  in  terms  of  vanishing  traction 
issues.  It  is  possible  to  obtain  a  robust  full  moment  tensor  solution  that  is  comprised 
dominantly  by  an  isotropic  or  explosive  component,  however  the  data  provide  the 
opportunity  to  evaluate  capabilities  of  moment  tensor  inversion  as  a  function  of  both 
source  depth  and  frequency.  The  details  of  the  data  processing  and  methodology  are 
similar  to  that  of  the  JVE  and  Lop  Nor  events,  as  described  in  section  3.2. 

Free-Surface  Vanishing  Traction 

The  first  step  before  performing  the  synthetic  analysis  is  to  look  at  how  the 
fundamental  Green’s  functions  behave  as  source  depth  decreases.  Using  the  western 
United  States  velocity  model  (Song  et  al.,  1996)  we  generate  the  ten  fundamental  Green’s 
functions  at  100-km  distance  with  source  depths  ranging  from  0.2  to  1.2  km.  The  Green’s 
functions  were  bandpass-filtered  between  10  to  50  seconds  period.  As  shown  in  Figure  13 
there  is  a  strong  source  depth  sensitivity  on  the  vertical  dip-slip  (DS)  fundamental 
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Green’s  functions  associated  with  the  Mxz  and  Myz  elements  for  all  three  components 
(ZDS,  RDS  and  TDS)  in  which  there  is  a  systematic  reduction  in  displacement  amplitude 
with  decreasing  source  depth.  This  effect  was  noted  in  a  study  on  fundamental  Love  and 
Rayleigh  waves  for  nuclear  explosions  and  associated  tectonic  release  (Given  and 
Mellan,  1986).  In  contrast,  the  vertical  strike-slip  Green’s  functions  for  all  three 
components  (ZSS,  RSS  and  TSS)  and  the  explosion  Green’s  functions  for  the  vertical  and 
radial  components  (ZEP,  REP)  show  little  to  no  variation  in  amplitude  and  waveform 
with  respect  to  depth.  The  45-degree  dip-slip  Green’s  functions  (ZDD  and  RDD)  show 
very  small  variations  in  waveforms  due  to  the  constructive  and  destructive  interference  of 
waves  interacting  with  the  free-surface.  While  the  wave  interference  appears  minor  in  the 
10  to  50  second  period  passband  (Figure  13)  it  is  more  pronounced  in  the  unfiltered 
displacement  Green’s  functions. 


ZDS  RDS  TDS 


Figure  13.  Green’s  functions  (GF)  computed  using  the  Song  et  al.,  (1996)  ID  western 
U.S.  velocity  model.  GF  are  in  displacement  and  filtered  between  10-50  second  period. 
The  traces  for  each  component,  from  top  to  bottom,  are  at  200,  400,  600,  800,  1000,  and 
1200  m  source  depths. 

The  low  amplitude  DS  component  Green’s  functions  could  lead  to  bias  in  seismic 
moment  tensor  results,  particularly  when  noise  in  the  data  is  also  considered.  It  is 
interesting  to  note  however  that  while  there  are  strong  effects  on  amplitude  the 
waveforms  remain  similar  and  there  is  little  effect  on  the  phase  of  the  waveforms  on 
these  components.  This  suggests  that  it  may  be  possible  to  develop  a  correctional  term  to 
scale  the  Green’s  functions  for  shallow  sources  prior  to  the  moment  tensor  inversion. 

In  the  synthetic  study  we  generated  a  suite  of  velocity  models  by  introducing  a 
shallow  velocity  gradient  to  the  ID  western  US  reference  model  (Song  et  al.,  1996).  This 
is  accomplished  by  splitting  the  top  2.5-km  thick  layer  in  the  reference  model  into  two 
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separate  layers.  We  systematically  adjust  the  thickness  and  velocity  of  the  two  new  layers 
(Figure  14),  but  constrain  the  variations  of  the  two  parameters  by  maintaining  the  same 
vertical  travel  time  as  the  reference  model.  The  purpose  of  the  study  is  to  generate 
different  but  comparable  velocity  models.  For  each  ID  model,  we  generated  Green’s 
functions  at  regional  distances  between  100  and  400  km  with  source  depths  ranging  from 
0.2  to  3.5  km.  Using  the  same  set  of  Green’s  functions  we  generated  two  types  of 
synthetic  data  with  different  source  mechanisms:  1)  a  pure  explosion  case  [EXP]  and  2)  a 
composite  source  case  (double-couple  [DC]  and  EXP)  with  an  F-factor  of  1  (Burger  et 
al.,  1986).  The  synthetic  data  are  computed  following  the  expressions  from  Minson  and 
Dreger  (2008).  We  added  random  Gaussian  white  noise  with  20%  of  the  maximum 
amplitude  of  the  synthetic  data,  then  we  use  a  10  to  50  second  causal  bandpass 
Butterworth  filter  to  filter  the  synthetic  data  and  the  Green’s  functions.  We  implemented 
a  semi-ideal  four-station  coverage  for  the  inversion,  consisting  of  source-to-receiver 
distances  at  100,  200,  300  and  400  km,  and  distributed  in  semi-regular  azimuths.  The 
linear  inversion  problem  yields  six  independent  components  of  Mij. 
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01234567  01234567 


Figure  14.  Velocity  models  derived  from  the  Song  et  al.  (1996)  ID  model  by  keeping  the 
top  2.5-km  vertical  travel  time  constant.  The  model  parameters  are  the  same  below  2.5- 
km  depth. 


Of  the  59  velocity  models  tested,  those  with  a  shallow  velocity  gradient  in  the 
upper  1.5-km  have  little  to  no  bias  in  their  isotropic  and  total  moment  estimates,  and  all 
estimates  fall  within  the  20%  noise  level  for  both  the  explosion  source  and  the  composite 
source.  Similarly,  moment  estimates  from  models  with  a  shallow  velocity  gradient  in  the 
upper  1.5-  to  2.0-km,  and  at  depths  greater  than  0.4  km,  are  all  within  20%  of  the  input 
values.  Although  the  total  moment  estimates  for  the  composite  case  exhibit  greater 
deviations  from  the  input  value  at  source  depths  shallower  than  0.4  km,  the  bias  in  the 
isotropic  moment  is  less  significant  (Figure  15).  In  most  cases  the  full  moment  tensor 
inversion  successfully  recovers  the  correct  mechanism  for  both  the  pure  explosion  case 
and  the  composite  case  over  the  targeted  depth  range  (<  1  km)  for  nuclear  explosions. 
Free-surface  vanishing  traction  has  little  effect  on  recovering  the  correct  mechanism  for 
models  with  a  shallow  velocity  gradient.  The  inversion  can  still  recover  the  correct 
mechanism  for  models  without  a  strong  shallow  velocity  gradient  at  depths  greater  than 
0.5  km;  the  bias  in  source  mechanism  is  significant  only  at  depths  shallower  than  0.4  km. 
The  key  observation  from  the  synthetic  tests  is  the  amount  of  bias  in  the  recovered 
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moment  and  source  mechanism  of  very  shallow  sources  depends  on  the  velocity  model. 
Models  with  a  shallow  gradient  produce  more  complex  Green’s  functions,  and  it  is 
possible  that  the  complexities  in  the  waveforms  provide  more  constraints  in  the  moment 
tensor  inversion  and  so  they  are  more  likely  to  recover  the  correct  solution. 


EXPLOSION  (EXP) 


Isotropic  Moment 


EXPLOSION  +  DOUBLE-COUPLE  (EXP+DC) 


Total  Moment  (Bowers  arid  Hudson,  1999) 


Depth  (km) 


CMT  Total  Moment  (Dziewonski  et  al.,  1981) 


Depth  (km) 


Depth  (km)  Depth  (km)  Depth  (km) 

Figure  15.  Isotropic  moment  and  total  seismic  moment  for  the  pure  explosion  source  and 
composite  source,  plotted  as  a  function  of  source  depth.  The  color  corresponds  to  the 
different  velocity  models  in  Figure  14.  Solid  black  lines  are  the  input/correct  seismic 
moments. 


HUMMING  ALBATROSS 

The  HUMMING  ALBATROSS  data  consist  of  both  broadband  and  short-period 
seismic  recordings  (Figure  16).  We  applied  the  moment  tensor  based  discrimination 
method  to  two  -M2-2.5  industrial  chemical  explosions.  Following  the  results  of  Ford  et 
al.  (2012)  and  Chiang  et  al.  (2014)  we  incorporated  full  waveform  data  from  five 
broadband  stations  and  P-wave  first  motions  from  sixteen  broadband  and  short-period 
stations  into  our  moment  tensor  analysis.  The  broadband  waveform  data  are  filtered 
between  0.5  to  2  seconds  depending  on  the  signal-to-noise-ratios,  and  we  examined  the 
results  in  two  frequency  bands  to  understand  the  frequency  dependence  of  free-surface 
vanishing  traction.  To  further  assess  the  uncertainties  in  the  moment  tensor  inversion 
result,  we  looked  at  how  the  solution  behaves  in  the  Hudson  et  al.  (1989)  source-type 
space.  Following  the  Network  Sensitivity  Solution  (NSS)  method  proposed  by  Ford  et  al. 
(2010),  we  randomly  generated  80  million  moment  tensor  solutions  that  are  uniformly 
distributed  on  the  Hudson  source-type  plot,  and  for  each  random  solution  we  compared 
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the  waveform  fits  and  P-wave  first  motion  picks  between  the  data  and  synthetics.  Here 
we  compare  two  types  of  NSS  result:  1)  waveform  data  only  and  2)  combined  full 
waveform  and  local  P-wave  polarity  data. 


O  Events 
■  Broadband 
A  Short  Period 


Figure  16.  Event  (shot)  and  station  locations  of  HUMMING  ALBATROSS.  Background 
colors  represent  topography  where  green  is  lower  elevation. 


Given  the  frequency  band  and  source-receiver  distance,  the  moment  tensor 
inversion  is  not  sensitive  to  source  depths  shallower  than  ~150  m,  but  the  solutions  at 
these  shallow  depths  remain  stable  and  are  dominated  by  a  large  isotropic  component 
(Figure  17).  This  is  true  for  both  shots.  If  we  have  no  prior  knowledge  on  source  depth, 
the  sensitivity  analysis  shows  we  can  constrain  the  depth  to  be  shallower  than  -150  to 
-200  m,  this  depth  range  is  indicative  of  possible  manmade  seismicity  since  natural 
earthquakes  rarely  occur  at  these  depths.  Because  we  know  the  depth  of  the  borehole 
where  the  explosions  are  detonated,  we  can  exclude  sources  deeper  than  15  m.  For  the 
subsequent  NSS  analysis  we  fixed  the  source  depth  at  1 1  m,  where  the  maximum  VR 
occurs,  even  though  the  VR  of  moment  tensor  solutions  between  5  and  15  m  are 
essentially  the  same. 
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Figure  17.  Full  moment  tensor  solutions  and  decompositions  as  a  function  source  depth. 
The  top  row  shows  results  from  shot  one  and  the  bottom  row  shows  results  from  shot  two. 
Note  the  x-axes  (depth)  scales  are  different. 


The  results  for  both  shots  are  very  similar.  Shot  one  has  a  slightly  larger 
magnitude  Mw  ranges  between  2.5  and  2.9  depending  on  the  frequency  band  used  in  the 
inversion.  For  Shot  two  Mw  is  between  1.9  and  2.  The  total  seismic  moment  increases 
when  we  use  longer  period  waveforms.  The  waveform  NSS  exhibits  a  tradeoff  between 
volumetric  sources  and  double-couple  sources  (Figure  18).  The  98%  VR  contours  wrap 
around  the  theoretical  DC  and  extend  towards  the  theoretical  opening  crack  mechanism, 
but  is  not  as  pronounced  for  shot  one  at  higher  frequencies  0. 8-2.0  Hz  (Figure  18a).  This 
tradeoff  is  likely  caused  by  the  combination  of  large  Love  waves  on  the  tangential 
component  and  free-surface  vanishing  traction.  The  large  Love  wave  causes  the  solution 
to  deviate  from  a  pure  explosive  mechanism.  Vertical  dip-slip  and  explosive  mechanisms 
fitting  equally  well  to  the  data  (Figure  19)  indicates  free-surface  effects  also  contribute  to 
the  observed  tradeoff.  The  tradeoff  becomes  more  pronounced  at  longer  periods  because 
free-surface  effects  are  more  severe  at  longer  wavelengths.  The  seismic  moment  also 
increases  when  we  look  at  longer  period  waves. 
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Figure  18.  Network  sensitivity  solutions  (NSS)  for  both  shots,  at  different  frequency 
bands  and  using  only  waveform  data  from  five  broadband  stations.  White  stars  are  the 
best-fitting  full  and  deviatoric  moment  tensor  solutions  from  time-domain  full  waveform 
inversion. 

Using  only  the  waveforms  we  cannot  constrain  the  two  shots  to  be  explosive  as 
the  uncertainties  in  the  NSS  encompass  non-isotropic  mechanisms.  To  eliminate  the 
tradeoff  we  applied  additional  constraints  from  P-wave  first  motion  polarities  (in  which 
all  show  upward  motion),  and  the  resulting  NSS’s  show  the  best  solutions  that  fit  both 
full  waveform  and  P-wave  first  motion  polarities  are  constrained  to  mechanisms  with  a 
dominant  isotropic  or  volume  increase  component,  including  the  full  moment  tensor 
solution  from  waveform  inversion  (Figures  19  &  20).  By  incorporating  additional  P-wave 
polarity  picks  from  the  local  network  of  stations,  we  now  have  greater  confidence  in 
discriminating  the  two  shots  as  explosive  events.  The  combined  NSS  results  are  equally 
well  constrained  for  both  frequency  bands.  Once  again  like  the  JVE  and  Lop  Nor  events, 
additional  P-wave  polarity  information  provides  a  very  important  constraint.  The  P-wave 
data  eliminates  mechanisms  such  as  the  vertical  dip-slip,  and  discriminates  the  event  as 
predominantly  explosive.  The  combined  results  show  that  the  moment  tensor  method  is 
capable  of  discriminating  very  shallow  events,  where  free-surface  effects  are  significant 
and  cannot  be  ignored. 
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Figure  19.  Combined  waveform  and  P-wave  first  motion  polarity  NSS  for  shot  one,  and 
the  full  and  deviatoric  moment  tensor  solutions  time-domain  waveform  inversion.  Red 
dots  are  upward  P-wave  first  motions  from  sixteen  broadband  and  short-period  stations. 
The  solid  black  lines  are  data  in  velocity  and  the  full  moment  tensor  synthetics  are  the 
dashed  red  lines. 
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Figure  20.  Combined  waveform  and  P-wave  first  motion  polarity  NSS  for  shot  two,  and 
the  full  and  deviatoric  moment  tensor  solutions  time-domain  waveform  inversion.  Red 
dots  are  upward  P-wave  first  motions  from  sixteen  broadband  and  short-period  stations. 
The  solid  black  lines  are  data  in  velocity  and  the  full  moment  tensor  synthetics  are  the 
dashed  red  lines. 
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3.4.  A  Continuous  Scanning  Algorithm  for  Micro-Earthquakes  and 
Shallow  Explosions 

We  have  applied  a  continuous  scanning  algorithm  to  determine  seismic  moment 
tensors  of  micro-earthquakes  with  no  a  priori  information.  This  method  may  be  applied  to 
streaming  data  in  real  time  to  detect,  locate  and  estimate  seismic  moment  tensors  of 
seismic  events.  It  may  also  be  applied  to  archived  data  to  investigate  large  numbers  of 
events  in  an  automated  manner.  In  the  following  we  describe  the  application  of  the 
algorithm  to  seismicity  associated  with  the  collapse  of  a  brine  cavern  on  the  side  of  a  salt 
dome  in  Louisiana.  These  events  are  unusual  and  display  large  volume  increase 
components  indicating  involvement  of  fluids  in  the  seismic  source  process.  The  results 
also  demonstrate  that  the  approach  could  be  used  to  detect,  locate,  determine  seismic 
moment  tensors,  and  discriminate  explosions  from  typical  earthquakes. 


Introduction 

Napoleonville  Salt  Dome  (. NSD )  is  located  at  -91.13°  E,  30.01°  N  near  Bayou  Come, 
Assumption  Parish,  southeast  Louisiana  (Figure  21).  It  is  a  part  of  the  Gulf  Coast  salt 
basin  which  exhibits  many  salt  structures  formed  by  upward  flow  of  sedimentary  salt 
(primarily,  evaporites)  on  account  of  low  density  of  salt  and  overburden  pressures  caused 
by  younger  sedimentary  deposits  (Beckman  and  Williamson,  1990).  It  penetrates  the 
Mississippi  River  Valley  Alluvial  Aquifer  ( MRAA )  zone,  which  is  primarily  composed  of 
upper  Pleistocene-Holocene  sediments,  to  an  approximate  depth  of  200  m  (Beckman  and 
Williamson,  1990).  Salt  domes  are  commonly  used  for  solution  mining  of  salt;  and 
caverns  thus  formed  are  also  used  for  storage  of  hydrocarbons  and  industrial  wastes 
owing  to  strong  impermeability  of  salt  (Thoms  and  Gehle,  2000).  54  caverns  distributed 
over  an  area  of  1.6  km  (N-S)  to  4.8  km  (E-W)  have  been  operating  in  NSD  at  various 
times  since  1950,  both  for  brine  mining  and  storage.  Beginning  in  June  2012,  residents  of 
Bayou  Come  reported  frequent  tremors  and  unusual  gas  bubbling  in  local  surface  water 
bodies  (Report  1,  2013).  The  parish  requested  the  assistance  of  the  United  States 
Geological  Survey  ( USGS)  to  monitor  the  continuous  seismic  activity.  A  temporary 
network  of  seismic  stations  was  set  up  which  revealed  a  sequence  of  numerous  seismic 
events  (Ellsworth  et  al.,  2012;  Horton  et  al.,  2013).  On  03  August  2012,  a  large  sinkhole 
was  discovered  close  to  the  western  edge  of  the  salt  dome  (Figure  21b)  leading  to  a 
declaration  of  emergency  and  evacuation  of  nearby  residents  (Report  2,  2013).  The 
sinkhole,  filled  with  a  slurry  of  water,  crude  oil  and  debris,  swallowed  cypress  trees,  and 
its  geometry  has  been  actively  changing  ever  since,  with  recent  estimates  of  surface  area 
being  greater  than  20,000  m2,  and  maximum  depth  varying  between  ~  30  and  100  m  over 
time  (Report  1,  2013).  Subsidence,  intermittent  seismicity  and  bubbling  of  natural  gas 
have  been  observed  in  the  affected  region.  Preliminary  investigations  suggest  that 
possible  collapse  of  the  lower  sidewall  of  a  plugged  and  abandoned  brine  cavern,  Oxy 
Geismar  3,  might  be  a  potential  cause  of  the  sinkhole  (Report  1,  2013;  Report  2,  2013). 
It  has  been  hypothesized  that  the  collapse  fractured  to  the  surface,  creating  a  disturbed 
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rock  zone  which  provides  a  pathway  for  formation  fluids,  natural  gas  and  crude  oil  from 
deeper  strata  that  are  now  accumulating  in  the  sinkhole  and  the  surrounding  subsurface 
(Report  1,  2013). 
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Figure  21.  (a)  Location  of  the  study  region  (black  star)  relative  to  the  state  of  Louisiana 
(southeastern  USA).  Google  Earth  Images  (dated  12  March  2013)  show-  (b)  the  study 
region,  indicated  by  the  rectangle  and  expanded  in  (c),  at  the  western  edge  of  NSD  (1000 
ft  and  10,000  ft  contours  indicated  by  white  lines;  William  Ellsworth,  personal  comm., 
2012);  (c)  Locations  of  the  5  USGS  broadband  stations  used  for  waveform  inversion 
(white  triangles),  approximate  location  of  Oxy  Geismar  3  cavern  (white  square)  and  an 
average  point  location  of  the  sinkhole  (white  balloon). 


Historical  seismicity  in  the  state  of  Louisiana  includes  a  widely  felt  magnitude  4.2 
earthquake  (peak  intensity  VI  on  MMI  scale)  on  19  October  1930  near  Napoleonville, 
Assumption  Parish  (Stover  and  Coffman,  1993),  whose  epicenter  might  have  been  close 
to  the  present  location  of  the  sinkhole.  However,  this  region  (-91.16°E  to  -91.13°E,  30  °N 
to  30.025  °N)  has  experienced  little  or  no  seismicity  in  recent  years,  with  no  earthquakes 
reported  in  the  National  Earthquake  Information  Center  ( NEIC)  -  Preliminary 
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Determination  of  Epicenters  (PDE)  Bulletin  between  January  1973  and  April  2012. 
Therefore,  the  synchronous  nature  of  seismicity  and  development  of  the  sinkhole 
suggests  that  the  two  phenomena  are  most  likely  related.  In  this  study,  we  investigate 
source  mechanisms  of  seismic  events  at  the  sinkhole,  represented  by  a  general  2nd  order 
point  source  centroid  seismic  moment  tensor  (MT).  We  implement  a  grid  search 
approach,  GRiD  MT  (Kawakatsu,  1998;  Tsuruoka  et  al.,  2009)  for  automatic  detection, 
location  and  MT  inversion  of  these  events  using  available  seismic  wave  velocity  models 
for  this  area.  We  show  results  for  the  time  period  of  one  day  just  prior  to  formation  of  the 
sinkhole.  We  check  the  stability  and  reliability  of  the  MT  solutions  and  interpret  them  in 
terms  of  possible  physical  processes. 

Data  and  Moment  Tensor  Inversion 

Three-component  waveforms  were  continuously  recorded  at  40  samples  per  second 
by  a  temporary  USGS  network  consisting  of  Trillium  Compact  broadband  seismometers 
and  a  Reftek  RT130  digitizer  (Figure  21c),  which  enables  the  study  of  this  seismic 
sequence  in  great  detail.  We  examine  events  starting  from  19:00  hours,  01  August  2012 
(just  after  station  LA09  became  operational).  As  signals  at  station  LA06  (not  shown  in 
Figure  21c,  but  located  SEE  of  sinkhole,  at  -91.12858  °E,  30.00771  °N)  are  very  weak 
and  have  higher  noise  levels,  we  haven’t  used  its  data  in  the  waveform  inversion.  Figure 
22a  shows  the  5-hour  (19:00  -  24:00  hours,  01  August  2012)  3-component  raw 
seismograms  at  station  LA02.  A  cursory  examination  shows  more  than  30  small  and  large 
events  spanning  one  order  in  peak  amplitude. 

Three  component  velocity  waveforms  at  station  LA08  for  two  of  these  events  are 
shown  in  Figure  23.  The  records  primarily  show  strong  surface  waves.  Our  analysis  of 
multiple  events  indicates  that  the  waveforms  are  quite  similar  to  each  other  indicating 
closely  spaced  hypocenters  and  a  repetitive  source  process.  For  most  events  the  duration 
of  strongest  ground  motions  is  limited  to  5-15  seconds. 

Variations  in  the  seismic  properties  of  the  subsurface  along  the  specific  paths  are 
evident  from  different  characteristics  of  the  broadband,  ~  0.1  to  15  Hz,  vertical 
component  velocity  waveforms  (Figure  24).  Stations  LA08,  LA01  and  LA02  are  on 
sedimentary  deposits  surrounding  the  salt  dome  and  show  waveforms  which  are  very 
similar  to  each  other,  but  quite  different  from  waveforms  at  stations  LA 09  and  LAO 3 
which  are  on  sediments  over  the  salt  dome.  This  suggests  that  the  salt  dome  structure 
affects  seismic  wave  propagation  over  these  relatively  short  source-receiver  distances. 
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Figure  22.  (a)  Raw  5-hour  record  on  01  August  2012  at  station  LA02.  The  arrow  at  the 
top  points  to  the  event  TE1  analyzed  in  detail  in  the  following  sections. 

Figure  22.  (b)  E-W  displacement  at  LA02  in  the  frequency  range  used  in  this  study  for  the 
same  time  period  as  (a).  Gray  lines  indicate  the  overall  background  signal  level  (±RMS) 
where  RMS  (Root  Mean  Square)  =  1. 76x1  O'5  cm. 
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Velocity  waveforms  for  two  different  events  at  station  LAO 8. 


Figure  24.  Vertical  component  velocity  records  for  one  event  ( TE1 )  at  all  five  stations. 
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The  time-independent,  2nd  order,  6-component  general  MT  (Jost  and  Herrmann, 
1989;  Minson  and  Dreger,  2008)  is  routinely  used  to  describe  the  source  mechanisms  of 
seismic  events.  The  scalar  seismic  moment  and  therefore  the  Mw  are  determined,  and  the 
MT  describes  the  mechanism  in  terms  of  body  force  equivalents,  i.e.,  double  couples  and 
linear  vector  dipoles,  assuming  a  point  source  in  space  and  time.  It  can  be  decomposed 
into  isotropic  (volume  change)  and  deviatoric  (zero  volume  change)  components. 
Typically  earthquakes  are  adequately  described  by  a  double-couple  ( DC)  mechanism  in 
which  one  of  the  eigenvalues  is  zero  and  the  other  two  are  equal  and  opposite  in  sign, 
although  small  non-double-couple  solutions  are  often  found  due  to  noise  in  the  data, 
unaccounted  path  effects  (Panning  et  al.,  2001),  or  possible  complications  in  the  source 
process  involving  non-planar  rupture  (Julian  et  al.,  1998).  The  MT  also  enables  the  study 
of  seismic  source  processes  for  non-earthquake  events  such  as  nuclear  explosions  (Ford 
et  al.,  2009a,  b)  and  mine  collapses  (Ford  et  al.,  2008)  that  involve  volume  changes  that 
result  in  significant  isotropic  components  (ISO)  in  their  MT  solutions.  The  MT  inversion 
methodology  is  described  in  detail  in  Jost  and  Herrmann  (1989),  Julian  et  al.  (1998),  and 
Minson  and  Dreger  (2008).  Requirements  of  a  standard  MT  inversion  are  a  hypocenter 
and  a  sufficiently  accurate  velocity  model.  The  inversion  then  solves  for  the  MT  and  the 
source  depth.  The  emergent  nature  of  first  arrivals  in  waveforms  of  the  seismic  events  at 
the  Louisiana  sinkhole  (i.e.,  in  Figures  23  and  24)  makes  accurate  picking  off-wave  and 
5-wave  arrivals  difficult.  Moreover,  given  that  the  distances  are  very  small,  a  small 
uncertainty  in  arrival  times  will  translate  to  a  large  relative  error  in  event  location.  The 
available  ID  velocity  models  (Figure  25)  of  the  study  region  (William  Ellsworth, 
personal  comm.,  2012)  show  large  differences  between  seismic  wave  velocities  at  the  salt 
dome  and  surrounding  sedimentary  structure,  which  is  in  turn  reflected  in  the  broadband 
seismic  waveforms  as  already  explained  in  Figure  24.  The  sediment  velocity  model 
shows  smoothly  increasing  velocities  with  gradients  decreasing  with  increasing  depth. 
The  velocity  model  over  the  salt  dome  consists  of  steep  velocity  gradients  associated  with 
sediments  over  the  salt  dome  up  to  a  depth  of  -250  m.  The  salt  dome  is  considered  to  be 
a  fast  half  space  with  Vp  =  4.5  km/s  and  Vs  =  2.5  km/s,  surrounded  and  overlain  by  slower 
sedimentary  layers.  The  3D  nature  of  the  subsurface  is  bound  to  influence  estimates  of 
the  source  location  and  MT  solution  at  shorter  periods,  however  considering  two  ID 
velocity  models  as  a  first  approximation  is  a  reasonable  approach  if  long-period  data  are 
used.  The  P-wave  and  .S-wave  first  arrival  times  computed  using  these  velocity  models 
and  the  ray-tracing  algorithm  of  Zhao  et  al.  (1992)  are  found  to  be  weakly  sensitive  to 
depth,  leading  to  large  uncertainties  in  depth  estimates.  Moreover,  the  large  number  of 
seismic  events,  more  than  100  on  some  days,  makes  a  standard  event-by-event  analysis  of 
location  and  MT  impractical. 
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Figure  25.  Velocity  and  density  models  used  in  this  study.  Gray  and  black  lines 
represent  salt  and  sediment  models,  respectively.  For  velocities,  thick  and  thin  lines 
represent  Vp  and  Vs,  respectively.  The  thin  dashed  black  lines  are  original  smooth 
velocity  models  with  gradients,  which  were  discretized  into  layers  and  used  to  compute 
GF. 


In  light  of  the  above-mentioned  difficulties,  we  employ  the  grid-search  approach, 
GRiD  MT,  of  Kawakatsu  (1998),  which  continuously  scans  the  seismic  wavefield  and 
performs  MT  inversions  of  relatively  low  frequency  waveforms  assuming  virtual  point 
sources  distributed  over  a  3D  grid.  For  a  given  time  window  of  data,  the  source  location 
and  MT  solution  which  give  the  best  Variance  Reduction  ( VR ),  a  measure  of  normalized 
goodness-of-fit  between  observed  and  synthetic  waveforms  (Minson  and  Dreger,  2008), 
is  inferred  to  be  the  true  seismic  source.  This  approach  has  been  used  for  detection  and 
source  characterization  of  offshore  earthquakes  in  Japan  using  streaming  long-period 
waveforms  (Tsuruoka  et  al.,  2009)  and  has  been  modified  for  a  tsunami  early  warning 
application  for  megathrust  earthquakes  (Guilhem  and  Dreger,  2011;  Guilhem  et  al., 
2013).  Based  on  preliminary  analyses,  we  construct  a  grid  extending  from  30.007  °N  to 
30.0136  °N,  -91.1452  °E  to  -91.138  °E  and  depth  0.02  km  to  1.77  km  with  grid  spacing 
0.0006°  and  50  m  in  horizontal  and  vertical  directions,  respectively.  This  grid  spans  the 
entire  present  volume  extent  of  the  probable  source  regions,  the  sinkhole  and  the  Oxy 
Geismar  3  cavern,  which  extends  in  depth  from  ~  1.0  to  1.7  km  (Report  1,  2013).  Three 
component  broadband  velocity  waveforms  are  first  corrected  for  instrument  response  and 
then  integrated  to  displacement.  At  sufficiently  low  frequencies  the  seismic  waveforms 
become  less  sensitive  to  small-scale  heterogeneities  in  earth  structure  making  simplified 
and  approximate  velocity  models  applicable,  and  the  seismic  source  can  be  treated  as  an 
effective  point  source  in  space  and  time.  We  have  tried  a  range  of  filters  with  different 
comer  frequencies  and  pole-orders  and  found  that  a  causal  4-pole  Butterworth  bandpass 
filter  with  comer  frequencies  at  0.1  Hz  and  0.2  Hz  greatly  simplifies  the  displacement 
waveforms  while  clearly  distinguishing  signals  of  larger  events  from  the  noise  floor 
(Figure  22b);  this  filter  is  subsequently  applied  to  both  observed  and  synthetic 
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waveforms.  It  is  also  important  to  note  that  there  is  a  narrowband  weakly  damped 
harmonic  signal  at  0.4  Hz  in  many  of  the  events  that  precludes  using  shorter  periods,  but 
which  may  be  indicative  of  a  triggered  acoustic  wave  phenomena  within  the  brine  fdled 
cavity,  or  resonance  possibly  induced  by  unsteady  fluid  flow  through  a  tensile  crack 
(Chouet,  1986;  Kumagai  et  al.,  2002).  Detailed  analysis  and  modeling  of  this  harmonic 
signal  is  beyond  the  scope  of  this  paper,  which  focuses  on  the  MT  analysis,  and  will  be 
investigated  in  future  work. 

For  waveform  modeling  in  subsequent  sections,  seismic  paths  to  stations  LA01,  LA02 
and  LA08  are  assumed  to  conform  to  the  sediment  velocity  model  and  seismic  paths  to 
stations  LA03  and  LA09  are  assumed  to  conform  to  the  salt  dome  velocity  model.  For  VP 
>1.5  km/s,  density  values  for  sedimentary  layers  are  computed  from  Vp  values  using 
Gardener’s  rule,  p  =  1.74  Vp0'25  (Gardener  et  al.,  1974).  For  VP<  1.5  km/s,  density  is  kept 
constant  at  1.93  g/cm  .  Salt  density  is  assumed  to  be  2.2  g/cm  .  At  such  small  distances 
(<  2.0  km),  which  are  smaller  or  comparable  to  seismic  wavelengths  being  used, 
waveforms  are  expected  to  be  weakly  affected  by  anelastic  attenuation.  Therefore, 
regional  QP  and  Qs  values  are  assumed  to  be  200  and  100,  respectively.  The  effects  of 
different  values  of  Q  on  synthetic  waveforms  were  analyzed  afterwards  and  displacement 
amplitudes  at  these  low  frequencies  were  found  to  be  insensitive  to  Qs  values  down  to  25, 
with  QP  =  2  x  Qs.  The  smoothly  varying  ID  sediment  and  salt  dome  velocity  models 
were  discretized  to  layered  models  and  Greens’  functions  ( GF)  were  computed  for  all 
grid  point  depths  and  all  distances  from  0.01  km  to  2.15  km  with  precision  of  0.01  km 
using  the  frequency-wavenumber  integration  software  FKRPROG  developed  by  C.K. 
Saikia  based  on  the  methods  of  Haskell  (1964),  Dunkin  (1965),  Watson  (1970),  Wang 
and  Herrmann  (1980)  and  Saikia  (1994).  GF  are  filtered  in  the  same  way  as  data  and 
decimated  to  0.5  seconds  sampling.  A  25  second  time  window  is  extracted  from  filtered 
data  and  then  decimated  to  0.5  seconds.  Assuming  each  grid  point  as  a  virtual  source,  a 
full  6-component  ( Mxx >  Myy,  Mzz,  Mxy,  Myz,  Mzx)  MT  inversion  is  performed  using 
expressions  from  Minson  and  Dreger  (2008).  The  VR  is  calculated  for  each  inversion  and 
is  used  to  assess  goodness-of-fit  and  to  identify  the  best  fitting  solution.  The  time  shift  to 
select  the  next  time  window  is  0.25,  0.5  or  1.0  seconds  depending  on  the  value  of  best  VR 
from  previous  time  window  (<  35%  or  between  35%  and  55%  or  >  55%,  respectively). 
This  ensures  that  the  grid  search  traverses  noise  windows  faster,  which  typically  return 
poor  values  of  VR,  while  maintaining  the  origin  time  resolution  at  0.25  seconds.  Coarse 
sampling  of  waveforms  and  adaptive  time  shifting  make  the  algorithm  computationally 
fast. 

GRiD  MT  is  applied  to  data  shown  in  Figure  22.  The  results,  in  the  form  of  the 
change  in  VR  with  time,  are  shown  in  Figure  26.  After  assuming  a  threshold  VR  of  70%, 
which  is  well  above  the  mean  background  VR,  we  find  23  events  during  this  period.  We 
have  also  detected  some  other  unusual  signals,  both  long  period  and  tremor-like,  which 
are  usually  restricted  to  the  closer  stations  {LAO 3,  LA08  and  LA09 )  and  might  be  related 
to  tilts  or  some  near-field  deformation  in  the  source  region,  unsteady  fluid  flow,  or  local 
noise  sources  (signals  not  correlated  across  multiple  stations).  Since  broadband 
characteristics  of  these  signals  are  completely  different  from  signals  of  the  discrete 
seismic  events  we  are  focusing  on,  they  are  not  included  in  this  study.  Details  of  the  MT 
solution  and  corresponding  waveform  fits  for  the  best  event  in  this  time  period,  event 
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TE1  (indicated  in  Figures  22  and  26)  are  shown  in  Figure  27a.  This  event  was  located  at 
grid  point  -91.1422  °E,  30.0112  °N,  depth  0.47  km  and  centroid  time  01  August  2012, 
20:52:39.00  hours.  The  solution  fits  the  data  very  well  at  VR  84.8%  and  can  explain  most 
of  the  strong  radial  and  vertical  components.  We  find  a  large  volume  increase  component 
(ISO  72%)  in  the  full  MT  solution.  Following  the  definition  of  Bowers  and  Hudson 
(1999),  the  scalar  Seismic  Moment  (Mo)  of  the  events  is  calculated  as,  Mo  =  Mj  +  MD, 
where  Mj  =  |(mi  +  m2  +  m3)/3|  is  the  isotropic  moment,  MD  =  max(|mj  -  (mi  +  m2  + 
m3)/3|)  is  the  deviatoric  moment,  and  mi,  m2  and  m3  are  the  eigenvalues  of  a  general 
moment  tensor  M.  Event  TE1  has  a  scalar  moment  thus  defined  of  2.4e+18  dyne. cm, 
corresponding  to  a  Mw  1.53.  The  spatial  distribution  of  VR  in  Figures  27b  and  27c  shows 
that  our  centroid  location  is  fairly  well  constrained,  and  located  east  of  the  sinkhole,  at 
the  edge  of  the  salt  dome.  Despite  the  long  seismic  wavelengths,  the  use  of  three- 
component  complete  waveforms  provides  both  arrival  time  and  azimuth  dependent 
polarity  information  for  various  phases,  thereby  strongly  constraining  the  locations.  We 
suspect  that  these  events  are  occurring  within  the  salt,  but  at  this  point,  we  are  unable  to 
precisely  put  the  event  in  the  salt  or  the  surrounding  sediments  due  to  the  coarse  grid 
spacing  (~  60  m),  the  large  seismic  wavelengths  of  data  used  in  the  inversion,  and  the 
uncertainties  in  3D  geometry  of  the  salt  dome  and  seismic  velocity  structure. 


Figure  26.  Plot  of  VR  with  time  for  data  shown  in  Figure  22.  Peaks  with  VR  >  70% 
(dashed  line)  are  considered  to  be  probable  seismic  events.  The  arrow  at  the  top  points  to 
event  TE1  indicated  in  Figure  22.  The  gray  line  indicates  the  mean  background  VR  (18.9 
%)  for  this  time  period. 
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Figure  27.  (a)  Observed  (solid  lines)  and  synthetic  (dashed  lines)  0.1 -0.2  Hz 
displacement  waveforms  and  full  MT  solution  for  event  TE1;  R  =  epicentral  distance,  az 
=  azimuth  (°),  Max  Amp  =  maximum  displacement  amplitude  at  a  station. 

Figure  27.  (b)  Grid-search  results  for  full  MT  solution  of  event  TE1  shown  in  (a).  Squares 
show  VR  at  grid  points  at  470  m  depth.  Solid  and  dotted  lines  are  1,000 ft  and  10,000 ft 
depth  contours  ofNSD.  Depth  sections  ofVR  across  profiles  A-B  and  C-D  (dashed  lines) 
through  the  best  fitting  centroid  location  are  shown  in  (c).  Shaded  polygon  (below  the 
grid)  =  approximate  surface  extent  of  the  sinkhole  in  July  2013;  black  triangles  =  station 
locations;  black  diamond  =  Oxy  Geismar  3  cavern. 

Figure  27.  (c)  Depth  sections  across  profdes  A-B  and  C-D  showing  smoothed  variations 
of  VR.  Black  star  is  the  best-fitting  centroid  hypocenter  of  event  TE1  and  the  thin  dashed 
line  is  the  line  of  intersection  of  the  two  sections. 

Figure  27.  (d)  Values  of  best  VR  at  various  grid  point  depths. 

Figure  28  shows  the  spatial  distribution  of  MT  solutions  represented  by  P- wave  first 
motion  mechanisms  for  all  23  events  detected  during  the  time  period  of  Figure  22. 
Reflecting  the  similarity  in  waveforms,  mechanisms  for  all  events  are  very  similar  to  each 
other.  The  events  are  approximately  concentrated  at  the  western  edge  of  the  salt  dome, 
very  close  to  the  sinkhole.  We  have  also  analyzed  the  data  of  02  August  2012  up  to  19:00 
hours.  We  observe  a  drop  in  seismicity  after  02  August  2012  corresponding  to  the  day 
when  the  sinkhole  was  first  discovered.  Overall,  the  mechanisms  and  locations  are  very 
similar  to  those  shown  in  Figure  28. 
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Figure  28.  Mechanisms  and  locations  of  the  events  detected  in  the  time  period  shown  in 
Figure  22.  Meaning  of  other  symbols  is  same  as  in  Figure  27b. 

The  MT  solutions  show  a  dominant  volume  increase  component,  which  is  quite  the 
opposite  to  what  one  would  expect  in  a  collapse  environment,  if  the  energy  release  were 
purely  due  to  gravity  driven  collapse  alone  (Ford  et  al.,  2008).  We  calculate  the  two  MT 
source-type  parameters,  e  and  k ,  for  all  events  and  plot  them  on  the  Hudson  source-type 
plot  (Hudson  et  al.,  1989)  shown  in  Figure  30.  The  horizontal  axis  plots  the  ratio  of  the 
deviatoric  eigenvalues,  e,  and  the  vertical  axis  plots  k ,  the  relative  isotropic  component. 
MT  solutions  for  the  Louisiana  sinkhole  seismic  sequence  plot  somewhere  close  to  tensile 
cracks  and  explosions,  quite  far  away  from  natural  double-couple  earthquakes  and 
expected  implosions  or  closing  cracks.  Decomposition  of  full  MT  solutions  of  these 
events  returns  6-32  %  DC  and  0-27  %  CLVD  components  in  addition  to  the  spherical 
tensile  source.  While  the  DC  components  are  very  small,  the  corresponding  fault  plane 
solutions  are  remarkably  similar  for  all  events,  indicating  that  the  small  deviatoric 
components  in  the  full  MT  solutions  are  not  random  and  are  possibly  due  systematic 
source  or  path  effects.  The  two  mean  fault  planes  (along  with  standard  deviations) 
consistent  with  the  DC  components  of  the  MT  solutions  are  (Strike  229°±19°,  Rake  - 
65°±22°,  Dip  47°±7°)  and  (Strike  18°±7°,  Rake  -112°±19°,  Dip  53°±9°).  Using  the 
hypocenters  from  the  full  MT  GRiD  MT  results  (Figure  27b  and  27c),  we  also  perform  a 
grid  search  of  source  parameters  for  a  shear-tensile  (crack+DC)  MT  solution  for  event 
TE1  (Minson  et  al.,  2007).  A  shear-tensile  source  mechanism  combines  tensile  opening 
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with  shear  slip  along  a  single  fault  plane  (Minson  et  al.,  2007;  Sileny  et  al.,  2009). 
Assuming  a  Poisson’s  ratio  (v)  of  0.39  for  a  source  in  sediments,  we  obtain  shear-tensile 
seismic  scalar  moments,  M0, crack  =  1.71e+18  dyne. cm  and  M0,dc=  3.7e+17  dyne. cm  on 
fault  planes  (Strike  225°,  Rake  -58°,  Dip  67°)  or  (Strike  23°,  Rake  -120°,  Dip  76°),  which 
fit  the  waveforms  well,  at  almost  the  same  VR  (~  84.8  %)  as  the  full  MT  solution.  A 
source  in  salt  (v  =  0.28)  yields  a  solution  that  fits  the  waveforms  poorly  at  VR  ~  78.4  %. 
This  is  also  reflected  in  the  Hudson  plot  in  Figure  29,  where  the  sinkhole  events  plot 
much  closer  to  a  theoretical  tensile  crack  in  sediments  (v  =  0.39)  than  to  one  in  salt  (v  = 
0.28).  Keeping  the  hypocenter  fixed,  a  pure  crack  in  sediments  (v  =  0.39)  also  yields 
good  waveform  fits  at  VR  =  83.1  %  whereas  a  pure  isotropic  explosion  and  a  pure  DC 
source  fit  the  waveforms  poorly  at  VR  38.5  %  and  60.2  %,  respectively. 


Figure  29.  Mw  and  source  depth  distribution  of  all  62  events  detected  with  VR  >  70%. 

To  assess  confidence  in  the  source  mechanisms  thus  obtained,  we  compute  the 
Network  Sensitivity  Solution  (NSS)  (Ford  et  al.,  2010)  for  event  TE1.  The  NSS  compares 
fits  between  observed  and  synthetic  waveforms  for  a  large  population  (3xl07)  of  source- 
types,  uniformly  distributed  moment  tensors  that  generates  a  distribution  of  goodness-of- 
fit  in  source-type  space  that  allows  one  to  identify  the  uniqueness  of  the  source-type,  and 
the  existence  of  possible  tradeoffs  as  is  common  in  nuclear  explosions  (Ford  et  al.,  2010). 
For  each  coordinate  on  the  Hudson  source-type  plot  (e,  k),  the  best  VR  value  from  all  MT 
solutions  corresponding  to  that  coordinate  is  selected,  and  normalized  with  respect  to  the 
maximum  VR  from  the  entire  space  (Figure  30a).  It  is  observed  that  MT  solutions  that 
produce  best  fits  (>  95%)  are  clustered  tightly  in  a  region  between  theoretical  explosions 
and  tensile  cracks,  quite  far  away  from  theoretical  deviatoric  mechanisms,  which  produce 
fits  only  up  to  70-75  %  of  the  best  possible  VR.  To  estimate  uncertainties  in  MT  solution 
of  event  TE1,  we  compute  10,000  full  MT  solutions  by  bootstrapping  residuals  from 
waveform  fits  for  the  best  fitting  full  MT  solution  at  the  same  location  and  depth  (Ford  et 
al.,  2009a,  b).  The  bootstrapped  residuals  are  first  filtered  using  the  same  filter  applied  to 
waveforms  and  are  then  rescaled  so  that  their  peak  amplitude  is  equal  to  the  peak 
amplitude  of  the  original  residuals.  The  99%  confidence  ellipse  of  distribution  of  these 
solutions  in  Hudson  space  is  shown  in  Figure  31a.  The  distributions  of  MT  elements  and 
DC,  CLVD  and  ISO  components  for  these  solutions  are  shown  in  Figure  31b.  The 
standard  deviations  of  normal-like  distributions  of  the  MT  elements  represent  the 
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uncertainties  in  the  MT  solution.  The  uncertainties  in  the  larger  MT  elements,  Mxx,  Myy 
and  MZz,  are  of  the  order  of  ~  2e+17  to  4e+17  dyne. cm  which  is  ~  9  to  30  %  of  their 
absolute  values  (1.3e+18  to  2.1e+18  dyne.cm).  While  the  uncertainties  in  the  smaller  MT 
elements,  Mxy,  Myz  and  Mzx,  are  large,  about  30%,  50%  and  1000%  of  the  absolute 
values,  respectively,  their  absolute  values  and  maximum  range  considering  2  standard 
deviations  are  smaller  than  the  absolute  values  of  Mxx,  Myy  and  MZz  by  an  order  of 
magnitude,  thereby  strongly  constraining  the  mean  ISO  component  to  ~  70%  with  a  small 
uncertainty  (±  4%).  The  uncertainty  in  Mxy,  Myz  and  Mzx  makes  discerning  between  a 
vertical  tensile  crack  vs  shear-tensile  failure  uncertain.  The  shape  of  the  99%  confidence 
ellipse  estimated  by  bootstrapping  waveform  residuals  closely  follows  the  shape  of  the 
98%  contour  of  best- fitting  source-types  in  the  NSS,  indicating  an  overall  consistency  in 
these  estimates. 


Figure  30.  Hudson  source-type  plot  showing  major  theoretical  seismic  source 
mechanisms  (black  crosses),  tensile  cracks  in  various  media  and  62  events  of  the 
Louisiana  sinkhole  seismic  sequence. 

For  shallow  events,  Ford  et  al.  (2009b;  2010;  2012)  observed  that  CLVD  sources  with 
vertical  axis  in  compression  provided  similar  quality  of  fits  to  explosion  waveforms  due 
to  their  mimicking  of  explosion  source  radiation  pattern  near  the  equator  of  the  focal 
sphere  which  is  the  region  usually  sampled  by  shallow  source  intermediate-period 
waveforms.  In  this  case  we  do  not  see  the  same  sort  of  tradeoff  in  the  NSS  since  the  near 
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vertical  tensile  crack  with  a  primarily  horizontal  major  vector  dipole  produces  strong 
azimuthal  variation  in  amplitudes,  which  is  not  produced  by  explosions  or  vertical  CL  VD 
sources. 
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Figure  31.  (a)  Network  Sensitivity  Solution  for  event  TE1  using  waveforms  only.  The 
white  polygon  is  the  99%  confidence  ellipse  of  the  distribution  of  MT  solutions  computed 
by  bootstrapping  residuals. 

Figure  31.  (b)  Distribution  of  DC,  CLVD  and  ISO  components  and  MT  elements  (units 
are  le+20  dyne.cm)  from  the  bootstrap  uncertainty  analysis  of  event  TE1.  The  dashed 
lines  represent  values  corresponding  to  the  best  fitting  full  MT  solution. 
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Results 


In  this  section,  we  have  implemented  a  procedure  for  independent  continuous 
detection,  location  and  MT  inversion  of  seismic  events  at  the  sinkhole  at  NSD,  Louisiana. 
The  computational  efficiency  and  simplicity  of  the  approach  makes  it  suitable  for  real¬ 
time  applications  or  analyzing  large  volumes  of  micro-seismic  data  in  reservoir  and  mine 
settings,  especially  if  events  are  numerous  or  one  is  not  confident  in  hypocenter  estimates 
based  on  travel-times  alone,  as  well  as  for  monitoring  for  possible  explosive  events  in 
targeted  regions.  The  method  has  the  advantage  that  given  a  reasonably  calibrated 
velocity  model,  and  use  of  relatively  long  period  waves  (actual  period  depends  on 
distances  involved),  that  it  is  autonomous  providing  independent  event  detections, 
locations,  source  parameter  estimation,  and  when  considering  full  moment  tensor 
implementations  source  type  discrimination  capability.  This  work  has  been  published  in 
the  Bulletin  of  the  Seismo logical  Society  of  America  (Nayak  and  Dreger,  2014). 
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3.5.  Full  MT  Analysis  of  the  Berkeley  Moment  Tensor  Catalog  and 
Demonstration  of  a  Statistical  Discrimination  Basis 

A  question  the  current  research  addresses  is  whether  source-type  populations  remain 
separated  in  other  regions  of  the  world,  as  is  the  case  for  NTS  and  western  US  seismicity. 
We  have  compiled  a  database  of  naturally  occurring,  and  induced  seismicity  to  examine 
the  distribution  of  non-double-couple  seismic  moment  tensors  in  a  Hudson  source-type 
diagram.  This  database  is  to  be  used  as  an  a  priori  constraint  on  the  complementary 
probability  distribution  of  an  explosion.  In  addition,  we  have  computed  full  moment 
tensor  solutions  for  877  events  in  the  Berkeley  Seismological  Laboratory  catalog  (Figure 
32a).  These  solutions  have  all  been  reviewed  and  are  of  the  highest  quality.  The  mean  is 
28=0.00544  and  k=-0. 03499,  essentially  a  double-couple,  and  the  data  is  normally 
distributed  to  first  order.  The  results  show  that  the  explosion  and  large  positive  CLVD 
area  of  the  source-type  diagram,  typical  for  explosions,  is  a  low  probability  region  for 
natural  seismicity.  We  have  evaluated  the  significance  of  improved  fit  of  full  moment 
tensor  solutions  compared  to  deviatoric  moment  tensor  solutions  for  the  Berkeley  catalog 
using  the  F-test.  These  results  indicate  that  only  20  of  the  events  have  statistical 
significance  above  90%,  and  that  the  vast  majority  of  events  (97.7%)  are  best 
characterized  as  deviatoric  without  a  volumetric  term.  Of  the  20  events  with  F-test 
statistical  significance  above  90%  only  7  have  statistical  significance  above  98%.  These 
events  occur  in  the  Long  Valley  Caldera  and  the  Geysers  geothermal  field  where  events 
with  large  positive  volumetric  components  have  been  well  constrained. 

In  Figures  33  and  34  we  compare  the  earthquake  and  explosion  data  populations  to  a 
P-value  to  test  for  evidence  of  explosion.  The  explosion  population  included  only  events 
from  the  Ford  et  al.  (2009)  study  of  NTS  nuclear  explosions.  The  earthquake  population 
is  the  BSL  catalog  shown  in  Figure  32.  For  each  population  a  multivariate  normal 
distribution  was  fit  to  the  data.  The  distributions  were  then  integrated  to  find  the  P-value 
for  the  earthquake  and  explosion  null  hypotheses.  A  composite  P-value,  which  is  the 
product  of  the  P-value  from  the  NTS  explosion  data  set,  and  the  complement  of  the  P- 
value  for  the  BSL  full  moment  tensor  earthquake  catalog  was  estimated.  The  shading  on 
Figures  33  and  34  show  the  composite  P-value.  This  composite  P-value  tests  for  evidence 
of  an  event  not  being  an  earthquake,  and  simultaneously  for  it  being  an  explosion.  As 
Figure  33  shows  only  a  few  earthquake  events  plot  in  the  region  of  relatively  high  P- 
value  for  explosion  evidence,  only  6  for  P-value>0.5.  Thus  6  out  of  877  (0.68%)  test 
potentially  high  as  an  explosion,  however  only  two  of  these  have  F-test  significance  that 
would  indicate  that  the  isotropic  component  is  resolved.  These  two  events  occurred  in  the 
Geysers  geothermal  field  where  events  with  isotropic  components  have  been  previously 
observed.  Considering  a  P-value  threshold  of  0.1  there  are  19  out  of  the  877  events  (2% 
of  the  events)  that  misidentify.  The  majority  of  these  have  F-statistics  that  would  indicate 
the  isotropic  component  of  the  moment  tensor  is  not  statistically  significant.  The  events 
in  this  region  with  F-statistics  above  95%  (red  pluses)  are  located  in  the  Geysers. 
Misidentification  of  earthquakes  is  a  matter  for  concern.  These  events  also  tend  to  have 
few  stations,  and  thus  the  isotropic  component  is  not  well  resolved,  and  it  is  likely  that  if 
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independent  first  motion  data  were  incorporated  that  the  false  isotropic  component  would 
be  eliminated.  Examining  misidentified  earthquakes  is  the  focus  of  ongoing  work. 

In  Figure  34  the  explosion  dataset  is  compared  to  the  composite  P-value  test.  The 
green  plusses  are  the  NTS  events  used  to  determine  the  source  type  distributions  for  the 
P-value  analysis,  and  the  red  plusses  are  test  events  that  include  the  1988  Soviet  JVE, 
1995  and  1996  Lop  Nor,  1998  India,  2006  and  2009  DPRK  nuclear  tests,  and  a  forth 
shallow  explosion.  All  of  the  events  discriminate  clearly  except  for  the  1996  Lop  Nor 
(section  3.2)  and  1998  India  tests.  Both  of  these  tests  have  large  F  factors  indicating 
reversal  of  Rayleigh  wave  polarity.  It  is  noted  however  that  while  these  two  events  fail  to 
discriminate  clearly  given  the  composite  test  they  are  also  inconsistent  with  both 
explosions  and  earthquakes.  Thus  even  though  explosion  discrimination  fails  for  these 
events  they  are  flagged  as  unusual  requiring  additional  analysis. 

The  example  provided  here  is  promising  in  terms  of  developing  a  statistical 
discrimination  basis.  It  makes  use  of  the  Hudson  et  al.  (1989)  source-type  diagram. 
Recently  Tape  and  Tape  (2012)  proposed  a  spherical  projection  of  source-type  that  is 
preferred  for  developing  a  statistical  based  discrimination  method.  Future  work  will  be 
implementing  their  projection. 


Figure  32.  BSL  full  moment  tensor  catalog  containing  877  events  from  1992-2012.  The 
symbols  are  color  coded  to  the  F-test  statistical  significance  comparing  a  5-degree  of 
freedom  deviatoric  moment  tensor  solution  to  the  6-degree  of  freedom  solution  that 
includes  the  volumetric  term.  The  ellipses  show  the  50%  and  90%  population  density 
curves. 
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Figure  33.  The  BSL  full  moment  tensor  catalog  (red  plusses  for  F-test  significance 
greater  than  95%)  is  compared  to  the  P-value  for  explosion  evidence. 


Figure  34.  NTS  explosions  (green  plusses),  and  test  explosions  (1988  Soviet  JVE,  1995 
Lop  Nor,  1996  Lop  Nor,  1998  India  2006  and  2009  DPRK  nuclear  tests,  and  a  shallow 
conventional  explosion;  red  plusses)  are  compared  to  the  composite  P-value  to  test  for 
explosion  evidence. 
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4.  RESULTS  AND  DISCUSSION 

The  results  of  this  study  indicate  that  a  regional  distance,  complete  long-period 
waveform  moment  tensor  approach  is  able  to  uniquely  discriminate  explosions  from 
earthquakes.  The  application  of  this  method  to  the  2009  DPRK  event,  section  3.1, 
indicates  that  the  theoretical  tradeoff  between  a  pure  explosion  mechanism  and  a  CLVD 
mechanism  with  the  major  compressional  dipole  oriented  vertically  can  be  eliminated 
using  P-wave  first  arrival  polarities  from  teleseismic  arrays  (Figure  2cd).  This  work  was 
published  in  Ford  et  al.  (2012).  Although  polarities  provide  substantial  additional 
constraint  as  demonstrated  there  can  be  difficulty  in  reliably  picking  first  P-wave  onset 
for  shallow  buried  sources,  and  therefore  using  the  teleseismic  array  beam  waveforms  is 
expected  to  provide  stronger  constraint.  In  fact,  as  Figure  2  shows  the  forward  fit  to  the 
array  beam  P-waveforms  is  quite  good,  and  it  should  be  possible  to  incorporate  these 
waveforms  in  the  inverse  procedure  for  determination  of  the  seismic  moment  tensor. 

The  results  of  the  study  of  the  1988  Soviet  JVE  test,  and  two  nuclear  explosions  at  the 
Lop  Nor  test  site  reveal  that  a  combined  long-period  waveform  and  first-motion  polarity 
moment  tensor  inversion  is  able  to  discriminate  these  explosions  from  earthquakes  under 
conditions  of  very  sparse  recording  geometry  (as  few  as  three  stations)  with  stations 
located  at  100s  to  1600  km  from  the  source  region.  The  incorporation  of  regional  distance 
P-wave  polarities,  at  the  same  stations  used  for  long-period  waveform  constraints,  or  if 
available  from  teleseismic  arrays,  greatly  improves  the  resolution  of  source-type  (Figures 
6  and  7)  even  for  cases  in  which  the  tectonic  release  may  be  strong  enough  to  reverse 
polarity  of  the  Rayleigh  wave. 

Synthetic  testing  indicates  that  the  vertical  dip-slip  Green’s  functions  associated  with 
the  Mxz  and  Myz  moment  tensor  components  are  affected  by  the  vanishing  traction  at 
the  ffee-surface  in  the  period  range  we  are  interested,  between  10  to  50  seconds  period. 
The  amplitudes  of  the  Green’s  functions  decrease  systematically,  but  the  waveforms  are 
similar  over  the  targeted  depth  range  for  nuclear  explosions  with  little  phase  distortion. 
Our  results  show  the  degree  in  which  free-surface  vanishing  traction  affects  the  moment 
tensor  solution  depends  strongly  on  the  velocity  model.  Velocity  models  with  a  shallow 
velocity  gradient  show  little  to  no  bias  in  the  isotropic  and  total  scalar  seismic  moment, 
for  both  pure  explosion  source  and  composite  source  mechanisms.  Similarly,  we  can 
retrieve  the  correct  mechanism  for  these  models  using  the  full  moment  tensor  inversion. 
One  possible  explanation  is  that  models  with  a  shallow  velocity  gradient  result  in  more 
complicated  waveforms  with  numerous  phases  more  completely  sampling  the  focal 
sphere,  thereby  providing  more  constraint  on  the  moment  tensor  inversion. 

Our  studies  show  that  including  P-wave  first  motions  in  addition  to  full  waveform 
data  can  eliminate  the  common  ISO-CLVD  tradeoff  (e.g.  Ford  et  al.,  2012;  Chiang  et  al., 
2014)  and  reduce  the  uncertainties  of  sparsely  recorded  underground  explosions  with 
strong  Love  waves  (Chiang  et  al.,  2014).  Here  we  demonstrated  that  although  we  cannot 
uniquely  characterize  the  HUMMING  ALBATROSS  events  as  predominantly  explosive 
using  only  waveform  data,  the  combined  waveform  and  first  motion  method  enables  the 
unique  discrimination  of  these  small  shallowly  buried  events.  The  combined  method  not 
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only  applies  to  larger  manmade  and  natural  seismic  events,  but  also  small  magnitude, 
very  shallow  explosive  sources  that  are  effectively  at  the  free  surface.  The  combined 
method  gives  a  well-constrained  NSS  even  as  we  go  towards  longer  periods,  where  the 
effects  of  ffee-surface  vanishing  traction  are  more  pronounced.  The  combination  of  both 
low  frequency  full  waveform  data  and  high  frequency  P-wave  polarities  greatly  enhances 
the  capabilities  of  the  moment  tensor  source-type  discrimination  method  in  cases  of 
sparse  station  coverage,  strong  Love  waves  and  free-surface  effects  due  to  very  shallow 
source  depths. 

A  continuous  waveform  scanning  algorithm  was  implemented  to  study  micro¬ 
earthquakes  associated  with  the  collapse  of  a  salt  cavern  in  Louisiana.  This  application 
demonstrated  the  recovery  of  stable  full  moment  tensor  solutions  that  show  these  micro¬ 
earthquakes  are  dominated  by  volume  increase  components  due  to  the  migration  of  high 
pressure  gas  and  fluids  for  the  subsurface.  This  example  shows  that  analysis  of  such 
events  can  be  carried  out  to  low  magnitudes  (Mw  0.8  to  2.1),  and  for  sources  at  very 
shallow  depth  450m.  In  this  analysis  stations  were  located  between  0.5  to  2  km.  While 
this  application  was  for  mining  induced  seismicity  it  demonstrates  that  it  is  possible  to 
autonomously  detect,  locate,  and  estimate  seismic  moment  tensor  solutions  using 
streaming  seismic  waveform  data.  Thus  the  approach  can  operate  in  a  realtime  processing 
stream.  It  also  has  an  advantage  in  that  if  a  particular  region  of  interested  is  well 
calibrated  with  3D  velocity  models  that  it  would  be  possible  to  incorporate  3D  Green’s 
functions  into  the  processing  stream  with  no  reduction  in  run  time  since  Green’s 
functions  are  stored  as  an  in-memory  database.  Using  the  principle  of  reciprocity  it  would 
be  possible  to  incorporate  path  specific  Green’s  functions  for  each  of  the  virtual  sources 
considered  in  a  given  virtual  source  volume.  This  approach  is  being  used  routinely  in 
Japan  to  monitor  offshore  seismicity  (Tsuruoka  et  al.,  2009),  and  we  are  currently 
implementing  the  method  for  seismic  monitoring  of  the  Mendocino  Triple  Junction 
region  to  be  used  as  a  tsunami  early  warning  algorithm  (Guilhem  and  Dreger,  2011),  and 
appears  to  be  very  well  suited  for  targeted  nuclear  explosion  monitoring. 

Using  a  catalog  of  earthquakes  for  northern  California,  and  the  NTS  nuclear 
explosion  results  of  Ford  et  al.  (2009a)  we  have  demonstrated  a  statistically  based 
discrimination  metric  that  utilizes  the  complement  of  the  P-value  for  earthquake  sources 
(the  hypothesis  that  events  are  not  earthquake  like),  and  the  P-value  for  the  NTS  nuclear 
explosions  (the  hypothesis  that  events  are  explosion  like).  The  result  of  the  composite  test 
was  applied  to  the  India,  JVE,  Lop  Nor  and  DPRK  nuclear  tests.  All  events  cleanly 
discriminate  except  for  India  and  the  1996  Lop  Nor  tests,  which  show  large  Love  wave 
excitation  and  evidence  for  reversed  polarity  Rayleigh  waves.  In  these  two  cases  the 
events  fall  in  a  region  that  would  indicate  a  suspicious  source  signature. 
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5.  CONCLUSIONS 


This  research  demonstrates  the  utility  of  a  regional  distance,  long-period  waveform 
moment  tensor  inversion  approach  for  nuclear  event  discrimination.  In  section  3.1  it  is 
shown  that  the  combination  of  long-period  regional  distance  waveforms  and  teleseismic 
P-wave  polarities  from  array  data  may  be  combined  in  a  manner  that  results  in  the  unique 
identification  of  the  2009  DPRK  nuclear  test  as  an  explosion.  In  section  3.2  we  further 
develop  this  approach  and  apply  it  to  very  sparse  monitoring  conditions  for  the  1988 
Soviet  JYE  and  1995/1996  Lop  Nor  nuclear  tests.  In  these  cases  it  was  found  that  the 
combination  of  long-period  regional  distance  waveforms  combined  with  P-wave  polarity 
information  from  the  same  stations  uniquely  discriminate  these  events  even  under  the 
extremely  poor  azimuthal  and  distance  coverage.  Additionally,  it  is  demonstrated  that  for 
the  1996  large  F  factor  test  in  which  Love  waves  are  larger  amplitude  than  long-period 
Rayleigh  waves,  and  where  the  tectonic  release  is  large  enough  to  cause  reversal  of 
polarity  of  the  Rayleigh  waves  the  combined  data  approach  coupled  with  the  NSS 
method  (Ford  et  al.,  2012)  identifies  the  event  as  unusual  compared  to  nearby 
earthquakes.  In  section  3.3  the  effect  of  shallow  depth  of  burial  and  the  consequent 
reduction  of  Green’s  function  amplitude  due  to  vanishing  traction  at  the  free  surface  does 
not  have  a  detrimental  effect  on  the  recovery  of  seismic  moment  tensors  for  shallow 
explosions,  however  this  is  velocity  model  dependent.  The  presence  of  shallow  structure 
that  produces  multiple  phases  that  more  thoroughly  sample  the  focal  sphere  appears  to 
diminish  the  effect  of  underestimation  of  explosion  scalar  moment  and  moment  tensor 
bias  that  can  be  observed  in  cases  with  simpler  velocity  structures  used  for  Green’s 
functions.  An  application  with  a  shallow  quarry  blast  shows  that  the  combination  of  long 
period  waveforms  and  first  motions  leads  to  unique  discrimination  of  the  event  as  an 
explosion.  In  section  3.4  we  demonstrate  the  application  of  a  continuous  scanning 
algorithm  that  enables  the  simultaneous  detection,  location  and  estimation  of  the  seismic 
moment  tensor  from  streaming  waveform  data.  The  application  in  this  case  was  for 
shallow  micro-earthquakes  associated  with  a  salt  dome  cavern  collapse,  however  it  shows 
that  the  method  may  be  applied  to  very  shallow  sources,  and  can  provide  a  complete 
monitoring  capability  for  targeted  regions  provided  velocity  models  and  Green’s 
functions  are  suitably  calibrated.  In  section  3.5  it  was  demonstrated  how  seismic  moment 
tensor  information  may  be  used  to  define  a  statistical  decision  basis  for  nuclear  event 
discrimination.  The  statistical  metric  used  is  the  combination  of  the  P  values  for  a 
calibration  population  of  nuclear  explosions  and  the  complement  of  the  P  value  for  a 
calibration  population  of  earthquakes.  The  combined  metric  was  shown  to  discriminate  a 
number  of  nuclear  explosions  not  used  in  the  calibration.  Although  the  1998  India  and 
1996  Lop  Nor  tests  fail  the  discrimination  they  are  identified  as  suspicious  being  outside 
of  the  distribution  of  naturally  occurring  earthquakes. 
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